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ABSTRACT 

Short, high intensity laser pulses induce nonlinear optical effects in the atmosphere 
that have the potential to make them propagate for long distances. Applications for 
long distance propagation of short pulses include active spectral remote sensing and 
laser lightning control. Much of the work in this field has been done with infrared 
pulses; however, it has been proposed that ultraviolet pulses have the advantage that 
longer pulse lengths can be used, thereby delivering more energy. Long pulse lengths 
lead to a simplified instantaneous model for the plasma response, which has been 
shown by Schwarz and Diels to admit steady state or oscillatory solutions corre¬ 
sponding to beam propagation. We have verified this model and have adjusted it to 
achieve closer agreement with numerical results. 

In this work we investigate the effects of transient behavior, and the stability 
of these solutions. Analysis of the modulational instability is done from the plane 
wave level to a full three dimensional model of the propagation. It is shown that 
both the transient behavior arising from the finite pulse length, and the modulational 
instability cause pulses to fragment over lengths on the scale of meters. We present 
results showing the growth of unstable modes in various propagation regimes. We 
discuss the pertinent length scales for ultraviolet pulses, as well as the impact of the 
instability and transient effects on theory and experiment. The results imply that 
continuous-wave models are very limited when used to predict dynamical properties 


of pulse propagation. 
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CHAPTER 1 
INTRODUCTION 

1.1 Motivation 

The study of short light pulses in the atmosphere is motivated by the variety 
of nonlinear processes that affect short pulses. Low intensity, spatially wide laser 
beams diffract, and are adversely affected by scattering and atmospheric turbulence. 
However, short pulses that are intense enough to induce several nonlinear optical 
effects in the air may be able to propagate for long distances, by changing the optical 
properties of the air. 

The primary nonlinear effect is self-focusing, which is caused by the intensity 
dependent refractive index induced by the optical Kerr effect [TJ. If the peak power 
exceeds a certain critical power, this self-focusing overcomes diffraction and causes 
the beam to focus tighter as it propagates, which in turn intensifies the nonlinear 
effects. At some point, the field becomes strong enough to ionize the molecules in the 
air, and higher order effects such as multi-photon absorption, plasma absorption, and 
plasma defocusing stop the collapse M- The propagation of a light pulse depends 
on the interaction between self-focusing and the various collapse arrest mechanisms. 
If there is a stable balance between these effects, the beam may become self-trapped 
and form a filament OT- If the balance is unstable, the pulse may split into multiple 
filaments in both space and time [5|. 

Tightly focused light filaments have the potential to pass through atmospheric tur- 
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bulcnce and scattering media such as clouds [6] with less distortion than a continuous- 
wave (CW) beam encounters. This opens up applications in remote sensing. As an 
ultrashort pulse compresses temporally, it can generate a broad spectrum, or “white 
light supercontinuum” mm- This might be used in future Lidar (Light Detection 
And Ranging) remote sensing techniques for spectral analysis of pollutants or mali¬ 
cious compounds in the atmosphere mm- The plasma channels generated by light 
filaments may also have the ability to enable laser induced lightning [T2l 13] . With the 
development of high peak power pulsed lasers, experiments are now being done with 
atmospheric propagation of ultrashort pulses 11 Emu]. Propagation distances as 
long as 2 km vertical in the atmosphere m and 200 m horizontal in air na have 
been claimed; however, questions remain as to what conditions are optimal for long 
distance propagation. 

1.2 Ultraviolet Filaments 

Much of the theoretical and experimental work on filaments has been done in the 
infrared, typically with wavelengths near 800 nm and pulse lengths on the order of 
hundreds of fs. However, some researchers are investigating the possibility of creating 
self-trapped pulses in the ultraviolet wavelength region, near 248 nm. UV filaments 
have been studied and contrasted with those in the IR [3j. In the ultraviolet, the 
critical power for self-focusing in air is smaller than the critical power in the infrared. 
It has been proposed that one may be able to scale the length of a propagating 
filament up in time, to longer pulses on the order of hundreds of picoseconds to 
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nanoseconds, instead of ultrashort pulses on the order of femtoseconds p~6l ITT] . This 
was proposed because it was observed that the balance of self-focusing and plasma 
defocusing appears to depend only on the intensity of the held |16j. Additionally, the 
reduced contribution from avalanche ionization in the ultraviolet extends the upper 
limit on the pulse duration up to tens of nanoseconds, four orders of magnitude larger 
than in the infrared H7|. Scaling up to longer pulses could increase the amount of 
energy delivered. 

This has been investigated by Schwarz and Diels, who have developed an analytical 
approximation to the problem of propagating longer pulses in the UV mi- The term 
“long pulse” depends on reference scale. For UV pulses, a long pulse is one that is 
approximated by a beam that is independent of time, because the scale of the pulse, 
on the order of nanoseconds, is much longer than the time scale of the nonlinear 
response of the air. In their semi-analytical approach, a Gaussian ansatz for the 
spatial beam profile converts the partial differential equation for beam propagation 
into an ordinary differential equation for the evolution of the beam waist size. Their 
steady state approximation forces the nonlinear response due to the plasma to act as 
an effective yA) nonlinearity. This approximation should be valid for pulses that are 
long enough for the generated plasma density to reach a steady state value. Their 
model predicts trapped solutions which oscillate in width, depending on the initial 
conditions. 

In Chapter [3j we numerically verify the Schwarz-Diels semi-analytical results, by 
using a split-step Fourier method of beam propagation. We have also been able to 
adjust some of their approximations, to make the model predictions quantitatively 
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agree better with numerical simulations. We then obtain idealized steady state CW 
beam solutions to the nonlinear propagation equation. However, this model, which 
contains no time dependence, cannot represent transient effects, such as how the beam 
solution is created in time, or whether the beam solution is stable. 

1.3 Experimental Status 

Many experiments have been done to study infrared filaments, often near 800 nm, 
and using short pulse lengths j6j 8, 9, 14, 15, 18-25]. Experiments range from labora¬ 
tory studies of filament formation in various media to atmospheric propagation tests. 
One project of note is the Teramobile, which is a mobile terawatt-femtosecond class 
laser operated by a joint French-German team mm. This laser generates 5 TW 
peak power pulses with 70 fs duration at 790 nm. It fits in a cargo container, and is 
being used to demonstrate atmospheric Lidar applications. 

Fewer experiments have been done to demonstrate UV filaments. The pulse 
lengths involved in current UV experiments, while typically longer than those con¬ 
sidered ultrashort in the IR, are shorter than those that will be discussed in this 
dissertation. For example, Schwarz et al. have performed experiments using a fre¬ 
quency tripled ThSapphire laser to obtain filaments at 248 nm. They created pulse 
widths between 600 fs and several ps, and observed filaments that propagated for 
12 m, much longer than the Rayleigh range of 12.5 cm for their pulses [16], [26] . 

Tzortzakis et al. performed experiments for both 450 fs and 5 ps pulses. They have 
observed ultraviolet filaments at 248 nm that propagate stably for 4 m. They noted 
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that their filaments are more stable than corresponding filaments in the infrared [4], 

An unanswered question is what experiments will tell us about longer ultraviolet 
pulses, such as lengths from hundreds of picoseconds to nanoseconds. It has been pro¬ 
posed that the results for ultraviolet filaments should scale with respect to increasing 
the pulse length [T6l [26] . resulting in long distance propagation of long duration UV 
filaments. However, experimental evidence for this is not yet available. The exam¬ 
ination of the theory for longer UV filaments, on the order of hundreds of ps, will 
occupy much of the following chapters. 

1.4 Filament Stability 

The semi-analytical solution for long pulses, developed by Schwarz and Diels, can 
be related to variational methods that have been developed for nonlinear optical prop¬ 
agation. These variational methods have been used P3 122123 to provide analytical 
or semi-analytical solutions to difficult problems, while retaining physical insight into 
the nature of the solutions. They are often used in conjunction with a model in 
which the temporal response of the plasma is assumed to be instantaneous. This is 
accomplished by either using a long pulse approximation ra. or, for shorter pulses, 
by examining the integrated contribution of the plasma at an arbitrary central time 
slice of the pulse [T5j [30, 31]. However, we will see that the problem of temporal sta¬ 
bility is critically important, and may be overlooked by using such approximations. 
We will investigate the behavior of the UV filament solutions by moving beyond the 


steady state model. 
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To investigate this question of stability, we develop a linear perturbation analysis 
that allows numerical simulation of oscillatory perturbations to a proposed zeroth 
order solution. We also use a three dimensional numerical code that includes the 
time domain and two transverse spatial dimensions. By first using periodic boundary 
conditions for the time domain, we isolate the fundamental instability from effects 
caused by the edges of the pulse. The results indicate that these pulses do suffer 
from a modulational instability. In Chapters [4] and [6j further analysis and simulation 
results detail the nature of the instability: its growth rate, frequency structure, and 
spatial form. 

We will also consider the relationship between the pulse length and the length 
scales of the various effects. Schwarz and Diels established that the pulse length must 
satisfy certain conditions in order to use their long pulse approximation HZj. If a 
pulse is too short, the steady state plasma approximation becomes invalid PH- If 
a pulse is too long, avalanche ionization begins to become important and must be 
added to the plasma model m- Since the growth rate of the instability depends 
on the frequency of the perturbation, in general the pulse length can play a role 
in stability. If a pulse is long enough to contain the strongest unstable periodic 
variations, it will be most strongly affected. In Chapter [?J we examine these length 
scales, and the experimental and theoretical implications of this instability and other 
effects for long ultraviolet pulses in air. We will also briefly examine the impact of 
stimulated Raman scattering [1], as it can produce instabilities, based on the time 
delayed nonlinear response [321 S3]. 
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CHAPTER 2 

REVIEW OF THEORY 

2.1 Propagation Equation 

This section briefly reviews the equations for propagation of a light field under the 
relevant nonlinear effects and approximations of interest. This will also help establish 
notational conventions and units. We will use MKS units throughout, unless stated 
otherwise. Start with the Maxwell wave equation in a nonmagnetic material with no 
free currents or charges pQ: 


V x V x E 


1 d' 2 E 
c 2 dt 2 


1 <9 2 P 

e 0 c 2 dt 2 ’ 


( 2 . 1 . 1 ) 


Assuming that V • E = 0, use the relation VxVxE = V(V-E) — V 2 E to obtain 


1 d 2 E _ 1 d 2 P 
c 2 dt 2 e 0 c 2 dt 2 ’ 


( 2 . 1 . 2 ) 


assuming a scalar field E and scalar polarization P. This reduces to 


n\ d 2 E _ 1 d 2 P NL 
c 2 dt 2 e 0 c 2 dt 2 


(2.1.3) 


where n 0 is the background linear index of refraction, if the material is assumed to 
be isotropic. In the quasi-monochromatic approximation, the real electric field E for 
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a wave propagating along the £-axis is then represented by a complex amplitude £\ 


E(x,y,z,t) = - £(x,y,z,t)e' {t ’ "*> + c.c. 


(2.1.4) 


Note that the exponential term is sometimes defined with the opposite sign, as in 
Schwarz and Diels na. With this particular complex representation, in the SI system 
of units, the field intensity is 




|£| 2 , 


(2.1.5) 


where Hq is the linear refractive index, eo = 8.85 x 10 12 F/m is the permittivity of 
free space, fio = 4nx 10 — ' H/m is the permeability of free space [1], and c 2 = (eo/ 4 ))” 1 - 


Some authors omit the factor of 1/2 from Equation (2.1.4). This alternate convention 


changes Equation (2.1.5) by a factor of four. 


Since we use scalar fields, we treat the nonlinear susceptibilities as scalars. A 
complete treatment for the general case of tensor susceptibilities and how they con¬ 
tribute to the nonlinear refractive indices can be found in the literature [1, 34]. The 
nonlinear polarization may be expressed as 


P NL = f(\£\ 2 )E = e 0 A X E. 


( 2 . 1 . 6 ) 


If the nonlinear response is instantaneous, Ay will be a simple function of |£| 2 . The 
quantity Ay may depend on time through the plasma density, which gives it a mem¬ 
ory based on the previous field history, which makes it non-instantaneous. The self- 
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focusing term (or 77 . 2 ) will also be non-instantaneous if the Raman effect is consid¬ 
ered; for now, let it be instantaneous. In this case, the simple nonlinear polarization 
term P NL for self-focusing can be written as 

P NL = e 0 x {3) \S\ 2 E. (2.1.7) 


The full nonlinear polarization will contain terms from self-focusing, plasma defocus- 
ing and absorption, and multi-photon ionization losses. 


Following Schwarz and Diels ra, insert Equation (|2.1.4[) for the electric held into 


the wave equation in Equation (2.1.3), to obtain 


y 2 1 _Kh ) c i(kz-wt) _ }_£!_ [Ay 

Vj - + a v 2 n 2 a+2 ce “ .2 f)+2 \ > 


< 9 2 \ 


1 3 2 


dz 2 c 2 dt 2 J 


c 2 dt 2 


( 2 . 1 . 8 ) 


where is the transverse Laplacian. Next perform the z and t derivatives above. 
Assume that £(x, y, z, t ) and Ay (a;, y, z, t ) are slowly varying in z and t compared to 
exp [i(kz — ut )\; neglect the second derivatives in z and t. This gives 


_q / d£ n 0 d£ . 

+ M (- + -—] + 


2 „, 2 


uj n \ 


0 


k 2 £ 


0 i{kz—uit) _ 


-2iw d ^* 8 ^ -uj 2 A x £ e i{kz ~ ut) . (2.1.9) 


Here we have used k = nouj/c and are neglecting any dispersion effects, such as group 
velocity dispersion (GVD). This allows cancellation of the last term in brackets in the 
above equation, and we group the d£/dz and d£/dt terms together. GVD and other 
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correction terms, which account for chromatic dispersion, are introduced through a 
series expansion [35] in studies of ultrashort pulses. However, we will either study 
the case of CW beams, in which there is no time dependence, or long pulses on the 
order of hundreds of picoseconds, for which neglecting dispersion is a very appropriate 
approximation due to their small spectral width. We will also neglect the d(A\E)/dt 
term on the right hand side. This correction to the nonlinear response is small for 
fields that are slowly varying compared to the optical frequency. Cancelling the 
exponential dependence gives 


For a CW beam, the time derivative of the held envelope is zero. However, if the 
held envelope is a pulse, we can use the transformation t' — t — (no/c)z, keeping the 
same z, to eliminate the time derivative. In both cases, rearranging terms gives 


d£_ 

dz 


i „ 9 „ ik . „ 

2k Vl£ + 2 ^, Ax 


( 2 . 1 . 11 ) 


For simple self-focusing, Ay = y^|£| 2 . Using the relationship y® = 2non 2 [H], this 
gives the self-focusing nonlinear Schrodinger equation 




( 2 . 1 . 12 ) 


where n 2 is the self-focusing index. Here y®, rt 2 , and \£\ 2 can all be in intensity or 
amplitude units, as long as one is consistent; see Appendix [A] for details. 
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If the nonlinear polarization includes the effect of the generated plasma, the prop¬ 
agation equation that we use is [16J 

+ ik 0 n 2 \£\ 2 £ - ^\£\ 2K ~ 2 £ « |(1 + iur)p£. (2.1.13) 

Here, the third term on the right hand side is a loss term corresponding to multi¬ 
photon ionization (MPI) generation, where j3^ is the coefficient for MPI at order K. 
The plasma term is derived from the Drude model [36]. In this term, a is the cross 
section for inverse bremsstrahlung, u is the optical reference frequency of the light, r 
is the characteristic electron collision time, and p is the electron plasma density. 

Note: In the paper by Schwarz and Diels H3. to which we often refer, there 
are some minor errors in units, discussed in Appendix [A] that do not affect their 
conclusions, but can cause confusion when comparing results. For simplicity, we use 
intensity units for all quantities such as £, and n^- Thus 77 - 2 / or 77 , 2 1£| 2 is the 
unitless quantity representing the refractive index change induced by the optical Kerr 
effect. 

2.2 Self-Focusing and Other Nonlinear Effects 

Some of the earliest work on self-focusing was motivated by the observation of thin 
damage streaks in materials in which an intense laser beam had been focused ra. 
The diameter of these damage tracks was smaller, by up to two orders of magnitude, 
than the corresponding diameter of a linearly focused Gaussian beam; the length of 
the tracks extended over several centimeters, indicating the presence of an intense 
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filament in which the light is tightly guided over many Rayleigh ranges 

The self-focusing is accounted for by the 112 term in the polarization, which de¬ 
scribes the optical Kerr effect [38]. This term is an intensity-dependent refractive 
index PQ. If n 2 is positive, a sufficiently strong held generates a region of higher index 
of refraction in the center of the beam, which acts as a focusing lens. The focusing 
further intensifies the held, which in turn increases the self-focusing effect p0]. For 
most of this dissertation, the 77,2 term will be assumed to describe an instantaneous 


response; in Section 7.4 we will discuss how including the Raman effect introduces a 
memory into this response, and the implications of this for filaments. 

When a beam contains more power than a certain critical power, it will undergo 
self-focusing collapse. As the beam intensifies, it begins to ionize the molecules in the 
air through multi-photon ionization (MPI). Collapse will continue until some physical 
effect, such as defocusing or losses from the plasma, acts to stop the collapse. The 
critical power for self-focusing of a Gaussian beam with wavelength A is 


p 

± r n 


A 2 


27m 0 n2 ’ 


( 2 . 2 . 1 ) 


where n 0 is the background index and 77,2 is in intensity units EH SO]. I 11 the ultraviolet 
at A = 248 nm, 77-2 is 7.8 x 10 -23 m 2 /W, and the critical power in air is approximately 
125 MW. In the near-infrared, the critical power is somewhat larger, due to smaller 
77-2 and larger A. At 775 nm, for example, the critical power in air is 1.7 GW [0U]. For 
pure self-focusing of a beam, the behavior only depends on the power in the beam, 
not its size or intensity. Beams with multiple critical powers present can break up 
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into multiple filaments due to modulational instability [5j SH 21 HIJ • 

Other relevant nonlinear effects appear as higher order terms, such as 77 . 3 1£| 3 , 
?r 4 |£| 4 , ..., which can act as defocusing terms if their sign is negative. However, 
terms that are odd in the electric held, such as the 77.3 term, only generally have 
nonzero coefficients in non-centrosymmetric media [TJ page 41], 

Intense light pulses generate plasma through multi-photon ionization of the mol¬ 
ecules in the air. This adds a nonlinear loss for the electric held. The plasma also in¬ 


duces nonlinear effects; this nonlinear response from the plasma, in Equation (2.1.13), 
includes both a plasma defocusing term and a plasma absorption term [13 Eg. For 
the special case of the long UV pulses discussed by Schwarz and Diels D3. the plasma 
defocusing was shown to act as an effective 77.3 term, even though air does not have an 
773 response. They also assume that the plasma response is instantaneous, compared 
to the time scales of long pulses. 

The balance between self-focusing, diffraction, and the other nonlinear effects has 
been the subject of many theoretical studies. One debate was whether the observa¬ 
tion of light filaments, which are extended in space, was due to the creation of an 
effective waveguide (self-guiding) or a the self-focusing of the various parts of a pulse 
(moving focus model). In the self-guiding model, a pulse is assumed to propagate 
stably under the influence of a wave guide created by the balance between the ion¬ 
ized plasma and the Kerr effect [2D] 02]. The moving focus model applies to finite 
duration pulses: it assumes that different time slices of a pulse will contain different 
peak powers, and thus focus at different distances, explaining the observation of an 
extended filament [24, 03]. In a third approach, the propagation of a short pulse 
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is described by a dynamic spatial replenishment model, where the pulse undergoes 
multiple collapses in time as it decays and is replenished by a new pulse pa mug. 

When we consider long ultraviolet filaments, we will first assume that they can 
be treated as sections of CW beams, and apply the long pulse model that will be 


described in Section [2T4j In this context self-guiding is the applicable model. 

More detailed discussions on the theory of self-focusing can be found in Boyd [Ij 
and Marburger [3lJ. Another, recent review of self-focusing in the context of wave 
collapse is provided by Berge 


2.3 Variational Techniques 


The equations for nonlinear beam or pulse propagation must in general be solved 


numerically, even when the approximations in Section 2.1 are made. An alternative 
approach is to use a variational method. One first casts the propagation equation in 
a Lagrangian form, and makes an ansatz for the form of a solution to the propagation 
equation. Then the Rayleigh-Ritz method or Kantorovich method [HI HI jHj, is 
used to determine the evolution of the parameters present in the ansatz. If the ansatz 
captures the physics of the solution well, it provides a simplified analytical or semi- 
analytical approach to solving a complicated problem. 

For example, Anderson and Bonnedal developed a variational treatment for non¬ 
linear self-focusing of Gaussian beams with a general nonlinear refractive index [28] . 
Their treatment includes a discussion of non-steady propagation, such as the oscil¬ 
latory solutions that Schwarz and Diels find. Wright et al. used their method to 
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study propagation of a CW beam in a medium with a cubic-quintic nonlinearity m 
Anderson and Bonnedal start with the equation 

d£ 

2ik— + + Q{\£\ 2 )£ = 0, (2.3.1) 

where Q represents a general nonlinearity. They present variational results for three 
forms of Q, the most applicable here being the series 

OO 

Q( |£| 2 ) = ]TC2,,|£| 2 “ (2.3.2) 

n =1 

By using the trial function 

£(r, z ) = Sq{z) exp[— r 2 /2a 2 (z) + ir 2 b(z)], (2.3.3) 


they obtain the variational result for the functional evolution of the beam radius (28j : 


d 2 a 
dz 2 


1 

k 2 a 3 


[1 - H(a 2 )] , 


where for the above form of Q, 


H(J) = E 

n= 1 


nC 2n (|go|ao) 2n 
(n + l) 2 a 2n ~ 2 


(2.3.4) 


(2.3.5) 


In this fashion, the partial differential equation for the evolution of the held £ is 
reduced to an ordinary differential equation for the evolution of the variational pa¬ 
rameter a. We will use their result later in Section |3f2| and compare it with the results 
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of Schwarz and Diels PH. Their method is not formally presented as a variational 
approach, but its results are nearly equivalent to other Gaussian beam variational 
approaches [27, ESI ;29J. 


We show in Section |3.2| that using the variational method for UV filaments in¬ 
dicates that the characteristics of the balance depend on the peak intensity, which 
agrees with existing treatments pm nz]. With just n 2 self-focusing present, the fo¬ 
cusing characteristics depend only on the power, for a Gaussian beam [33]. However, 
when the balance between n 2 and n 3 is analyzed, we find that the peak intensity of 
the beam is a more suitable indicator of its characteristics. For trapped solutions, 
the peak intensity, once chosen, determines the width and thus the peak power of the 
steady state solution, as well as its stability characteristics. This intensity dependence 
motivated the proposition that scaling to higher pulse energies and longer pulses in 
the UV might be possible mm 

Other authors have applied variational methods to pulse propagation, with in¬ 
creasing levels of detail. Akozbek et al. have studied femtosecond pulse propagation 
using a variational method similar to that discussed above [29] - They include the 
nonlinear loss terms through a dissipative modification of the standard variational 
treatment mm- 


2.4 Long Pulse Model 

In the long pulse model, or near steady state approximation, used by Schwarz and 
Diels, the action of the plasma on the light field is reduced to a simple instantaneous 
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nonlinearity. During the middle section of the pulse, the time dependence of the 
plasma is removed, by assuming it has had enough time to reach a steady state value. 
This will be a good approximation if the plasma generation rise time is much shorter 
than the time scales of the pulses of interest. The plasma term in the propagation 
equation is then replaced by an effective nonlinear index term. In the case of UV 
filaments at 248 nm, the plasma acts as an effective 713 term. The advantage of such 
a model is that if it is accurate, one no longer needs to consider the time dimension 
when doing a simulation, which becomes difficult for long pulses. 

However, as applied to the UV problem, the long pulse model also places an upper 
limit on the length of the pulse. The upper limit is set by inverse bremsstrahlung, 
which causes avalanche ionization, and in the Schwarz-Diels analysis comes out to 
4-60 ns. The lower limit, as mentioned, is set by the time it takes for the plasma to 
build up to its steady state value, and is is roughly 30-200 ps PHi26j. These limits 
depend on the strength of the multi-photon absorption rate, which is a value that is 
known with very poor accuracy, within only two or three orders of magnitude. As the 
MPI rate decreases, the upper limit decreases and the lower limit increases, narrowing 
the width of the region of pulse lengths that satisfy the model. 

Schwarz and Diels derive their semi-analytical solution for propagation of pulses 
with lengths that satisfy these limits. In the introduction to their paper, however, 
they mention that their results do not address the stability of the solution nor the 
transient evolution to steady state ra. In the following chapters, these issues will be 
shown to present a critical weakness of using a long pulse model for UV filaments. 
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CHAPTER 3 

STEADY STATE SOLUTIONS 

3.1 Application of the Long Pulse Model 

In this chapter, we review and extend the method of Schwarz and Diels S3, 
and cite their primary results to provide a foundation for the investigation into the 
stability and transient effects. Start with the following equation for propagation [16] 
of the electric held envelope S(x, y, z, t ) : 

zp = + ik ° n2 l^l 2 ^ ~ ^~^\ S \ 2K ~ 2S ~ \+ iuj r)p£. (3.1.1) 

The electron density p of the plasma is described by 

I = ^lD + ^|p-V + C v i A (3.1.2) 


where C is a coefficient for avalanche ionization, a is an electron-positive ion recom¬ 
bination coefficient m, and D is a diffusion strength. In steady state, the left hand 


side of Equation (3.1.2) is zero. If there is no diffusion (D = 0), the steady state 
plasma density is the solution of a quadratic equation at each point in space. 

Note that Schwarz and Diels mi use a different form for the multi-photon ioniza¬ 
tion source term in their equation for the plasma. Their MPI source term is a^I 3 No, 


where a^ is the three-photon MPI cross section, / is the intensity, and N 0 is the 
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density of oxygen in the air. Using the values in Table 3.1, note that 


^ = b = 1.6 X 10-“ DV 


3 hu 


W 3 s 


(3.1.3) 


This consistency is achieved by using the values in Ref. ra for cr® and Nq to derive 
the value of rather than using the value from Ref. [16], We will use the shorthand 
notation that b represents the overall MPI coefficient for the plasma. 

In the model used by Schwarz and Diels, the steady state ansatz allows one to 
express the plasma density p as a power of the held strength or intensity. This is 
only true for certain pulse lengths HZ]. First, if the pulses are short enough, less 


than tens of nanoseconds, the avalanche term in Equation (3.1.2) can be neglected. 


If the pulses are long enough so that the plasma has time to reach steady state, 
they set dp/dt = 0. This rise time is on the order of 30-200 ps, depending on the 
publication mm- They also neglect the diffusion term. Thus, solving for p, 


P = \fbja\£f, 


(3.1.4) 


where b = P^/Kfrw. Given the light frequency of interest, which tells us the MPI 
order K, under this approximation the plasma density is proportional to \E\ K . Thus 
the plasma response term can be replaced by an effective uk\E\ k £ nonlinear index 
term. For UV pulses at 248 nm, K = 3, and p oc \£\\ so Schwarz and Diels replace 
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the plasma term by an effective n 3 term. They derive the value of this term as 


n 3 = - 


a® N 0 


a 2noU} 2 m e eo ’ 


(3.1.5) 


where e and m e are the charge and mass of the electron, using a model for the electron 
plasma na. To reach the simplest case, neglect losses in the propagation equation, 


and consider them later once the solution is understood. Equation (3.1.1) reduces to 


d£_ 

dz 


+ ik 0 n 2 \£\ 2 £ + ik 0 n 3 \£\ 3 £. 


(3.1.6) 


Note that the plasma defocusing coefficient aurp/2 from Equation (3.1.1) is equal 
to kon 3 \£\ 3 , given the steady state dependence of p on |£| 3 , for the parameters in 
Table 13.11 


Schwarz and Diels use a semi-analytical approach with a Gaussian beam ansatz to 
convert this equation to one for an imaginary particle representing the beam width, 
moving on a potential surface. Their result is an ordinary differential equation that 
can be solved with a straightforward integration. If the particle is trapped in a 
potential well, the beam will oscillate in width as it propagates. They illustrate these 
oscillatory solutions in their paper HU; the size of the oscillations depends on how 
steep or shallow the well in the effective potential is. 

There is a question as to how the beam width oscillations are to be interpreted. 
Schwarz and Diels suggest mi that they correspond to a steady state analog of the 
dynamic replenishment phenomenon [40]. However, this model explains the propa¬ 
gation of short pulses in air by time-dependent interactions between the leading and 
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Tabic 3.1: Parameters used in propagation 


Parameter 

Name 

Value 

Units 

Background index 

no 

1 

— 

Self focusing index 

n 2 

7.8 x 10- 23 

m 2 /W 

Plasma nonlinear index 

n 3 

-3.3 x 10- 31 

m 3 /W 3 / 2 

MPI order 

K 

3 

— 

Wavelength 

A 

248 

nm 

MPI coefficient 


3.9 x 10 -34 

m 3 /W 2 

Recombination loss 

a 

1.1 x 10" 12 

m 3 /s 

Inverse bremsstrahlung cross section 

a 

5.2 x nr 25 

m 2 

Electron collision time 

T 

3.5 x 10- 13 

s 

Density of oxygen 

N 0 

5.4 x 10 24 

m 

3-photon absorption cross section 


3 x 10" 41 

m 6 s 2 /J 3 

Linear absorption loss 

OIL 

2.5 x nr 4 

m 1 


trailing edges of the pulse and the generated plasma; the pulse splits in time, and 
power from the trailing edge replenishes the pulse. The oscillations that we consider 
are for a time independent beam that can be characterized by a single width as a 
function of time. Additionally, any spatial pulse shaping, such as ring formation, 
cannot be described by the Gaussian beam ansatz. 

The parameters that we use in simulations are presented in Tabic |3.1| Many of 
these are from Schwarz and Diels dl, so that we can compare and verify the results. 
Additional parameters are taken from Ref. [16]. Some of the values in the table are 
not well known, in particular the MPI coefficient which determines the MPI rate. 
This MPI coefficient is difficult to measure [3Ej, and may change by two orders of 
magnitude up or down, depending on the publication [T6j [TT, 26] HE] • The numbers 
from Schwarz and Diels m, are values that will allow comparison to previous work, 
as well as give physically interesting insights into the application of the steady state 
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model. The following sections discuss details of the oscillatory solutions, and how 
the approximation in the Schwarz-Diels approach can be modified to give better 
agreement with numerical simulation. 

3.2 Improving the Aberrationless Approximation 

In Schwarz and Diels mi. the electric held amplitude £ (r, z) is assumed to be a 
Gaussian. Then, they make the following parabolic “aberrationless” approximation 
for powers of the electric held HU Equation 30], 



This approximation is substituted into the |£| 2 and |£| 3 terms in the nonlinear po¬ 
larization, and terms are matched based on their order in r 2 /w 2 . One consequence 
of this Taylor series approximation is that the critical power obtained using their 
analysis is one-fourth the critical power that is normally obtained in theories of self- 
focusing [271 HQ]. The factor of four comes from the difference between the true 
Gaussian beam and the parabolic approximation. 

However, other authors make a different aberrationless approximation, based on 
minimizing the mean square error 

(ci + c 2 r 2 — \£\ 2 ) 2 \£\ 2 2nr dr (3.2.2) 



in the polynomial approximation Ci + C 2 r 2 of the exponential |49j . With their method, 
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the aberrationless approximation for \£\ 2 becomes [49] 


-2 r 2 /w 2 



(3.2.3) 


This approximation will change the critical power to the accepted value. This form 
of the aberrationless approximation was used by Cerullo et al. to study space-time 
coupling effects in the collapse of femtosecond pulses [SD]. The weighting function for 
the mean square error should remain |£| 2 , independent of which power of \£\ is being 
approximated. When this is done for the plasma term, which depends on |£| 3 , the 
approximation becomes 


/16 12 r 2 \ 

\25 ~~ 25 w 2 / ' 


(3.2.4) 


The main results of the Schwarz-Diels analysis are the equations for the evolution 
of the beam width. We present these results here, with the above modification to 
their approximation. If we apply the new approximations to their analysis, beginning 
with Equation 32 of their paper El. their Equation 37, which we will show describes 
the potential surface, becomes 


dw\ 4 
dz ) k 2 


P 

W ~ 1 


ur 


Wr I 


— —/?4 

25 n 0 


■ur 


Wn 


p 3/2 + S. 

K 0 


(3.2.5) 


and their Equation 38, which is the differential equation for the evolution of the held 
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width, becomes 


£=4 (£ - A -4+V 2 

a,P k z \ P;, r J w 6 25 n 0 wr 


(3.2.6) 


where P c ' r = X 2 / (Zimonz) is four times the P cr in Schwarz and Diels |1?]. This gives a 
more accurate estimate of the true critical power for a Gaussian beam, which is very 
close to P' cr but not exactly equal mmm- The quantity /3 4 is defined in Schwarz 
and Diels as 


- 2 \ 3 / 2 

I'X - - ( J n 3 , 


(3.2.7) 


where n 3 = —3.3x 10 -31 m 2 /W 3//2 is the nonlinear index contribution from the plasma, 
in intensity units, and l = 8.9 mm is the mean free path length of the electrons in 
the plasma S3- 

This ordinary differential equation can be solved numerically, given the input 
beam width and curvature. Schwarz and Diels present results detailing the oscillatory 


solutions to this equation m In Section |3.3[ we will compare the solutions of this 
modified equation to results obtained in a numerical propagation of the beam. 

However, we can arrive at this same result with a different approach. Wright et 
al. [27] applied the general variational result of Anderson and Bonnedal [28], discussed 


in Section 2.3, to the problem of beam propagation under cubic-quintic (712 and 
nX) nonlinearity. This variational model assumes that the beam is described by the 


trial function in Equation (2.3.3). Following their notation, we can apply the same 


result to the Schwarz-Diels case of a cubic-quartic nonlinearity. In Anderson and 
Bonnedal [28], the nonlinear index terms only contain even powers of \E\ 2 . However, 
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their variational analysis is still valid for odd powers, by letting their series index 
take half-integer values. By considering the n$\£| 3 term, we find that the equation of 
motion for the beam width in the case of the UV filaments is 


d 2 a _ 1 / _ Pin 24 fTp P in a 0 \ 
dz 2 k 2 a 3 y P cr 25 V A P cr a J 


Here a is the beam width, and a 0 the initial beam width. Note that there are two 
different definitions of the beam width; w = \/2a, where w is the width that Schwarz 
and Diels use. P* n is the input beam power; P cr is the critical power as defined in 


Equation (2.2.1). I p is the peak input intensity, and J 0 is a reference intensity defined 


by Jq = (t^Ad) 2 - At I = Jq, the total nonlinear refractive index change for a plane 


wave is zero. 


For the numbers in Schwarz and Diels, Jo = 5.5 x 10 16 W/m 2 . Following Ref. 
we obtain a self-trapping power for a collimated Gaussian input beam 


Psr = Pa 


, 24/4 
25 V h 


-l 


(3.2.9) 


This self-trapping power depends on the input intensity. If I p > (25/24) 2 Jo, then 
there will be no self-trapped solution, since power cannot be negative. For UV 
propagation, this predicts that there will be no self-trapped solution for I p > 6.0 x 
10 16 W/m 2 . The self-trapped power, however, can be arbitrarily large, as the beam 
width may increase as intensity approaches the cutoff. 

If we convert the equation of motion from a to w and to the notation that Schwarz 
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and Diels use, we find 

<Pw = _± (Pin _ \ 1 24 

dz 2 k 2 \P cr ) w 3 25 n 0 w 4 


(3.2.10) 


Here we must note that P cr is the more accurate critical power, which is four times 
larger than that referenced in Schwarz and Diels. This result is exactly identical 


to that in Equation (3.2.6), which was obtained by modifying the aberrationless 
approximation, although it was obtained from a Gaussian beam variational ansatz. 


3.3 Verification of Oscillatory Solutions 

Schwarz and Diels used an analytical approximate approach to convert the partial 
differential equation for field propagation into an ordinary differential equation for 
an effective particle in a potential. They then numerically integrated the ordinary 
differential equation. We directly solve the numerical partial differential equation for 
propagation, using a split-step Fourier method, which is detailed in Appendix |B} 

In Ref. [IT], the field is described by its Gaussian beam width, which oscillates in 
an effective potential well. The potential well model for nonlinear beam propagation is 
often used and is a direct consequence of using a variational technique [27113111151116]. 
To understand the analogy of the potential well, recognize that 

KS) 2+f/w=s (3 - 3 ' i) 


is the equation for the position w of the particle with total energy E as a function of 
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the coordinate z, which acts as a classical time, as it evolves in the potential U{w). 


Equation (3.2.5) can then be written in the form 




(3.3.2) 


where 


TT ( n 2 ( P ^1,8 Z/3 4 P 3/2 
^ ° ) w 2 + 25 n 0 w 3 


fc 2 VP. 


is the effective potential and 


(3.3.3) 


E = U(wq) 


Wq_ 

2 Rl 


(3.3.4) 


is the total energy, which is conserved. The initial width wq of the beam gives the 
particle’s initial location; the initial width of the beam divided by its initial radius of 
curvature gives its initial velocity. This can also be seen from Equation 37 of Schwarz 


and Diels na, or Equation ( |3.2.5[ ) of this dissertation, by setting w = wq. Thus the 
initial velocity corresponds to any initial linear focusing or defocusing with which the 
beam is prepared. 

If the initial conditions are near the bottom of the well, the oscillations in the 
width are small. However, if the initial conditions displace the particle far from the 
well minimum, there are larger oscillations. The potential U(w) approaches zero as 
w —» oo. If the particle has total energy E greater than zero, it is not bound and 
the beam will diffract. This can be caused by the initial beam being too narrow, or 
having too much initial curvature. 
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We use the Schwarz-Diels semi-analytical solution with improved aberrationless 
approximation described earlier, to achieve closer agreement between the direct nu¬ 
merical simulation and their predictions. As an example, Figure [37T| shows the poten¬ 
tial surfaces corresponding to beam powers of 100 MW, 200 MW, and 500 MW. The 



Beam width (microns) 


Figure 3.1: Potential wells for beam powers at 100, 200, and 500 MW. 


surface for 200 MW is shallower than the surface for 500 MW. The 200 MW surface 
represents states which are more weakly bound. The 100 MW surface is for a power 
below the critical power of 125 MW, so the potential surface no longer supports a 
bound state. The location of the bound state minimum depends on the power as well, 


according to Equation (3.3.5). Since the binding strength of the potential depends 


on the power in the beam, absorption losses will cause the potential to change, and 
we will observe this in the oscillatory solution. 


To solve the modified evolution described in Equation (3.2.5), we use the ODE 
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solver in Mathematica. The initial condition for the radius of curvature of the beam 
is dw/dz|. =0 = w 0 /R 0 . For simplicity, start with a flat beam (R 0 = oo), so the initial 
width is also the waist. The location of the equilibrium, or steady state, beam waist 


in Equation (3.2.5), which is the modified version of Equation 37 of Schwarz and 


Diels, depends on power as 


w eq (P) = 


0.24 k 2 l(5 A P 3 ' 2 

n 0 PI Per- 1' 


(3.3.5) 


In the case that we wish to model absorption, the evolution of the beam power, in 
Equation 18 of Schwarz and Diels, is unchanged by the modification, so we reproduce 
it here: [T7] 

1 rlP P 2 P 3 / 2 

(3.3.6) 


1 dP _ a P 2 a P3 ' 2 

~ — P3^ ~ P4“ 


P dz r '°w 4 ^ w 3 


During the split-step numerical simulation, we characterize the beam width by a 


second moment computation, described in Equation (5.5.1). An example of the case 


where the initial conditions are near the bottom of the well is shown in Figure 3.2 


Here the power was 500 MW and the initial width 100 /iin, with no absorption. 
The frequency of the simulated oscillations roughly matches that of the predicted 
oscillations, but the prediction does not account for the observed decay in amplitude. 
For other values, such as those shown in Figure |3.3 the initial conditions lead to 
larger oscillations. In the second case the initial width was chosen to be 200 /iin. In 


the field history plot, Figure 3.3(b), the beam grows in intensity and decreases its 


width, while shedding rings during the oscillations. Figure |3H| shows the evolution of 
a beam where the power is 200 MW and initial width is 100 //m. Here there are fairly 
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consistent, larger amplitude oscillations, but again they are not constant. For starting 
powers larger than 1 GW, the beam quickly sheds much of its power into rings, so 
that the peak intensity is reduced. This is consistent with the existence of a maximum 
trapping intensity, although the observed value for this maximum trapping intensity 
is roughly 4 x 10 16 W/m 2 , instead of 6 x 10 16 W/m 2 as predicted in Section 
The variational approach provides an estimate for the order of magnitude of such a 
quantity, but cannot capture the details of features such as rings. 

The oscillations appear to either dampen in amplitude or change chaotically for 
these cases. The predictions, however, are based on a potential model with no loss 
terms, in which the particle should oscillate with a constant amplitude between two 
extremes. The apparent loss in the potential model should not be confused with 
nonlinear propagation losses, which we have not yet considered. The damping is 
caused by the fact that the evolving solution has to adjust itself, possibly by shedding 
power into rings, to fit the nonlinear response and match a solution of the propagation 
equation. If the initial condition is far from the equilibrium solution, the Gaussian 
beam approximation will be temporarily violated as the beam changes form. Any 
power that is diffracted in a ring will be absorbed by absorbing boundary conditions 
at the edges of the grid. After some propagation, the damping reduces, as the evolving 
solution changes from an initial Gaussian to a solution that oscillates in an stable 
fashion. The oscillatory propagation in these results is similar to that observed in 
studies of a medium with cubic-quintic nonlinearity [27, 45]. 

If propagation losses are present, the power in the beam is reduced, which leads to 
a shallower potential. If the power loss is not extreme, it should act as an adiabatic 


3.2 
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(a) Beam width using split-step evolution 



Distance (m) 

(b) Predicted beam oscillations using modified Schwarz-Diels model 

Figure 3.2: Beam evolution for power of 500 MW, initial width of 100 /irn, 
and no absorption. This set of initial conditions corresponds to a near 
steady state solution. 
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(a) Beam width using split-step evolution, together with regular 
oscillations predicted by Schwarz-Diels model 



(b) Early portion of field history using split-step evo¬ 
lution. The plot shows a cross section of the field 
amplitude, as a function of propagation distance. 


Figure 3.3: Beam evolution for power of 500 MW, initial width of 200 /irn, 
and no absorption. In this case the initial width is far from steady state, 
leading to nonlinear oscillations. 
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(b) Predicted beam oscillations using modified Schwarz-Diels model 


Figure 3.4: Beam evolution for power of 200 MW, initial width of 100 /nrl, 
and no absorption. 
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change in the potential surface. The oscillations should become larger in amplitude, 


but lower in frequency, as power is lost. This is shown for one case in Figure 3.5 


For this simulation, we included MPI losses and plasma absorption using the values 


from Table |3.1| The beam power appears to decay roughly exponentially, but since 
the loss is nonlinear, the loss rate decreases as less power remains. As power is lost, 
the effective potential becomes shallower; one might expect the oscillations to thus 
increase in amplitude. However, they do not do this; they are instead somewhat 
chaotic. They will disappear once the beam power drops below the critical power; 
then the potential is unbound, and the beam width continues to increase as it diffracts. 
The particular nature of the oscillations in this figure may be influenced by numerical 
effects such as reflections from the absorbing boundary conditions. However, the 
important point is that the beam width follows a general trend dictated by the slowly 
decreasing power in the beam. 

Recall that starting with a 500 MW Gaussian beam that is 100 /irn wide, for 
the ri '2 and n-j parameters in Schwarz and Diels ra, we find that the beam is close 
to an exact steady state solution. The location of this initial condition is near the 
bottom of the well on the potential surface for 500 MW shown in Figure [3H] However, 
even choosing an initial condition with the beam width at the precise bottom of the 
potential well will not yield us an exact steady solution, because the potential picture 
is approximate, based on a Gaussian ansatz; the exact solution is not a Gaussian. 
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Distance (m) 

(a) Beam width using split-step evolution 



(b) Beam power 


Figure 3.5: Beam evolution for 500 MW initial power, initial width 
of 100 /mi, and absorption losses. The beam width oscillates around a 
changing equilibrium position. The loss rate decreases for longer distances, 
dne to the nonlinear nature of the absorption. 
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3.4 Construction of Exact Steady State Solution 

For completeness, and for diagnostics that will be a part of the stability analysis 
discussed later, we wish to construct exact steady state trapped beam solutions to 
the propagation equation, within numerical precision limitations. Let a steady state 
solution be denoted by £o(r, z ). 

To accurately find £ 0 , we investigate the solution of the eigenvalue equation that 
corresponds to the steady state solution. The exact solutions will be characterized 
by their peak intensity. By solving this boundary value problem with the shooting 
method, we obtain a radial profile of £ 0 . This radial profile is then extended to a 
full x-y grid using interpolation. We look for a steady state solution to the simplified 
propagation equation 


^ = 2k v ' 2±£ + ik ^\ £ \ 2 + n *\ s ?) s - (3.4.1) 

Let 

£{r,z) = £{r)e ifiz . (3.4.2) 

with £{r) real. Then, 

i(3£{r)e ll3z = — V]_£(r)e^ z + ik 0 [ n 2 £(r ) 2 + n 3 S(r ) 3 ] £(^r)e %l3z , (3.4.3) 


and 


V\£{r) = 2 k/3£(r) - ik 0 [ n 2 £(r ) 2 + n 3 S(r) 3 ] £(r) = f(£). 


(3.4.4) 
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Recall that we are using intensity units for £, so £ 2 is an intensity. In polar coordinates 


d 2 £ \_d£_ 
dr 2 r dr 


m- 


(3.4.5) 


For a given peak amplitude Aq, the solution will have the initial conditions that 
£(r = 0) = Aq, and that dS/dr| r=0 = 0. Then, given a guess is for /?, we solve the 
ordinary differential equation, starting at r = 0 and running until r is several times 
larger than the extent of the expected solution. 

This equation caused difficulty for the numerical integration routine in Mathe- 
matica, clue to the singularity in the equation at r = 0. However, a finite difference 
approach can be used to handle the singularity. We represent the value of £(r) by 
values on a discrete uniform grid a 0 , cq,..., a at. The a m are values of £(r) at the 
points r = mA, where A is the grid spacing in r. 

The operator is approximated by the following difference formula: 



d £ 1 d£ &n +1 2a n T tt n — l 1 Un+1 1 

dr 2 r dr A 2 nA 2A 

n 

(2 n + l)a n+ i — 4na n + (2n — l)a n _i 
“ 2nA 2 ' 


(3.4.6) 

(3.4.7) 
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and we then derive an relation which can be iterated: 


4n 2n — 1 2 n 2 

° n+1 = 2^+i an - 2^+i (ln - 1 + 2 ^+ 1 A 


(3.4.9) 


This iteration equation is valid for n > 1. The initial conditions present a slight 
difficulty. We know that ao = 7L 0 , but we need to know a value for a±. The second 
derivative at n — 0 must be handled with care. To represent V]_, which is a 2D 
operator, on a 2D spatial (XY) grid using a finite difference technique, the value for 
V 2 is approximated by 


V 2 £ = 


Gffio — 2 oq,o + 0 - 1,0 a o,i — 2ao,o + ®o,-i 


(3.4.10) 


where now the two indices for a represent the X and Y directions on the grid. We do 
not define this function on a full XY grid, but we can approximate this difference as 


V 2 £ = 


4a i — 4a 0 


o A 2 


(3.4.11) 


since, due to the azimuthal symmetry, the four nearest neighbors all have the value 
a\. Then we can get the condition for ap 


Vi4 = = f( A o), 


(3.4.12) 


A 2 

a,i = a 0 + —f (A 0 ). 


(3.4.13) 
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With the initial values for ao and cl\ now set, we guess an initial value for f3 and then 
use the “shooting” method to adjust f3 until the value of the solution £{r) for large 
r approaches zero. In practice, once the guess for [5 is good, the solution will decay 
to zero at some distance, but for sufficiently large r, the solution then diverges from 
zero due to lack of precision. However, if the distance for which the solution diverges 
is larger than the grid we need in propagation, we can truncate the solution earlier. 
This approach gives the solution for S{r ) on a uniform grid in r. Then the solution 
is interpolated from the radial grid to a 2D grid to generate £(x,y). This solution is 
saved, and will be used as the steady state held £ 0 in the stability analysis and three 
dimensional propagation in later chapters. This exact steady state beam solution for 


a peak intensity of 3.2 x 10 16 W/m 2 is shown in Figure 3.6 We chose this value for 
the intensity because 3.2 x 10 16 W/m 2 is also the approximate peak intensity for a 
Gaussian beam with power 500 MW. The value of /3 for this solution is approximately 
8.11 m" 1 . 

We first test the generated So solution by propagating it in the 2D propagation 
code with the correct conditions for and no- The code has two ways of handling 
the plasma term, either by using an effective no, or the term — acorp/2. Recall that 
we chose the values for the the effective no and the rates used to generate p (the 
/ 3^ value) to be consistent. When the correct solution for £ 0 is loaded, the field 
stays constant in the propagation code, up to one part in 10 5 over a distance scale of 
several meters, if no losses are present. Thus the steady state solution also provides 
a verification of the propagation code, as well as the two different approaches to 
handling the plasma term in the code. 
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Figure 3.6: Example of numerically constructed steady state solution. 
This solution contains a power of roughly 470 MW, and the beam width 
is roughly 97 /mi. The peak intensity, which defines this solution, is 
3.2 x 10 16 W/m 2 in this example. Also shown is a Gaussian profile, with 
the same peak intensity and width. 


The steady state solution for this power is visually close to a Gaussian, but not an 
exact match. For pure self-focusing, the self-similar collapse solution is the Townes 
profile mmm-, in that case the difference between the Gaussian and Townes 
profiles is important, due to the unstable nature of the solution. For the UV problem 
under the Schwarz-Diels model, an initial Gaussian may be nearly steady state if 
chosen correctly. This indicates that the Gaussian ansatz in the Schwarz-Diels model 
can be a good guess, when describing steady state or near steady state solutions. As 
we saw, it is not as appropriate when describing non-equilibrium solutions. 

The results from this section confirm the oscillatory solutions predicted by Schwarz 
and Diels, and demonstrate that the simplified propagation equation has exact steady 
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state solutions for CW beams, up to a phase change as the beam propagates. 

The next question is whether long ultraviolet pulses can be sent through air, 
maintaining their shape over long propagation distances, and whether the current long 
pulse model can describe this. To answer this question, we must further scrutinize 
the steady state model representation of long laser pulses. As mentioned earlier, 
Schwarz and Diels impose upper and lower limits on the pulse length to make the 
approximation valid, so first the pulse must satisfy these constraints. The steady 
state model also assumes that the plasma responds instantaneously to changes in the 
electric held. Thus it has no way to incorporate any temporal transient effects or 
instabilities that may arise from the finite response time of the plasma to the electric 
held, and this is mentioned in the analysis im In the following chapters, we develop 
a stability model, and simulate the propagation of long pulses to investigate these 


effects. 
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CHAPTER 4 
STABILITY 

4.1 Introduction 

The stability of solutions to nonlinear optical propagation equations has been 
studied for many cases. The question of stability has been considered for the mathe¬ 
matical balance in a beam solution between self-focusing and diffraction, and this is 
known to be an unstable situation [DE5U32I21E3]. Only the mathematically exact 
solution with precisely one critical power will propagate unchanged; any deviation 
of the beam power from the critical power results in the beam solution collapsing or 
diffracting indefinitely. 

More detailed models consider the spatiotemporal instability of the solutions to the 
nonlinear Schrodinger equation (NLSE) [51 1541 [55] . Liou et al. studied the spatiotem¬ 
poral instability of an NLSE that only contained GVD and self-focusing terms [ 53] , 
They found that under self-focusing conditions, the interaction of spatial and tempo¬ 
ral modulations led to an instability for either normal or anomalous GVD, which led 
to spatiotemporal chaos developing in the solution. 

The addition of some form of saturation, such as the generation of plasma, stops 
filaments from collapsing indefinitely. The balance between the self-focusing collapse 
and saturation has been thought to explain the long life of the observed light filaments, 
or “bullets” that are created in beams with power equal to several critical powers [56]. 

However, the delayed responses of the Kerr effect and the plasma defocusing have 
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been shown to create instabilities in short (single picosecond scale) pulses [32] . Kan¬ 
didov et ah attributed the hlamentation or breakup of short pulses to these instabil¬ 
ities 


Bian et al. studied the instability of plasmas generated by laser pulses 153, which 
they denote the ionization scattering instability. They identified a dispersion relation 


for the growth of unstable modes. Our analysis in Sections 4.3 and 4.4 is similar to 


their approach; however, they consider the impact of this instability for short (less 
than a picosecond) pulses. 

Since many of these investigations are done for short time scales, the instabilities 
are identified as the breakup mechanism for a short pulse into filaments [21 HEMES]. 
However, this raises the question of the stability of the UV long pulse solutions, which 
are based on an instantaneous Kerr and instantaneous plasma defocusing response. 
Additionally, the Schwarz-Diels model assumes that the leading edge of the pulse 
can develop into the steady state middle region. This requires that the steady state 
solution be stable, since the leading edge, which deviates from steady state, has 
to eventually adjust itself into a steady state solution for the model to be valid. 
Our analysis will first use a linear perturbation model to investigate this question of 
stability. 


4.2 Purely Spatial Perturbations 

In general, we would like to understand the growth of unstable modes as a func¬ 
tion of both spatial frequency k± and temporal frequency, which we denote by 0. 
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First, examine the instability for zero temporal frequency, that is, a purely spatial 
modulational instability for a CW solution. We use an approach similar to that of 
Cerda [33]. Consider the propagation equation 

Pjf i 

f^ = ^VlS + ik 0 A n (\S\ 2 )E, ( 4 . 2 . 1 ) 

where An(\S | 2 ) is a general instantaneous nonlinear index change. For the long pulse 
UV problem 

An(\£\ 2 ) = n 2 \£\ 2 + n 3 (\£\ 2 ) 3 / 2 . ( 4 . 2 . 2 ) 

An exact plane wave solution is 

So(x,y } z) = S 0 e ikoAnz , ( 4 . 2 . 3 ) 

where An = An(|£o| 2 )- For the cubic-quartic nonlinearity in the UV problem, the 
propagation constant, denoted k, for the exact plane wave solution is then 

k = k 0 An = k 0 n 2 \£ 0 \ 2 + k 0 n 3 \£ 0 \ 3 . ( 4 . 2 . 4 ) 

For a beam solution, this plane wave value provides a magnitude estimate, but not 
an exact prediction, for the value of the propagation constant k. 

Consider a perturbed solution 


£(x,y,z) = £c,e lkoAnz (1 + € + (x,j /,z) +e.(x,y,z)), 


( 4 . 2 . 5 ) 
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where e + and e_ are small perturbation fields. Note that to first order 


|£| 2 ~ |£ 0 | 2 + |^o| 2 (e+ + €- + c.c), (4.2.6) 


so to first order, we can approximate 


An(|£| 2 ) « An(|^ 0 | 2 ) + |£o| 2 (e+ + e_ + c.c) 


dAn 

d\£^' 


(4.2.7) 


Insert the expression for £ into Equation (4.2.1) for the propagation. After collecting 


terms, we obtain equations for the evolution of the perturbations: 


de + - i Y7'2 ,-T. dAn \C 12/ , m 

2k v±t+ + lk °a\sp ] o{ (e+ ° 

<9eA f 2 9An 2 * N 

“ _ 2*: v±£+ “ |£o1 (e+ + '->• 


(4.2.8) 


Assume that the perturbation fields are plane waves such that 


e + (x, y, z) = u + exp(Az + • r) 

e_(x, y, z ) = «_ exp(A*£ — ?kj_ • r). 


(4.2.9) 


This ansatz for the form of the fields, when inserted into the coupled evolution equa¬ 
tions above, gives a matrix equation: 


—T + V V 


u+ 

= A 

u + 

l 

A 

1 

h 

A 

1 


u*_ 


u*_ 


(4.2.10) 
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where 

r) /\ 7i 

T = k 2 j2k, V = k 0 —\£ 0 \ 2 . (4.2.11) 

The eigenvalues A satisfy 

A 2 = 2TV-T 2 = T(2V-T). (4.2.12) 

Note that T and V are both real, and T > 0. If V < 0, then A 2 < 0, so A is 
imaginary, and the system is stable. But if V > 0, A 2 > 0, so it is unstable for 
T < 2V. T corresponds to the transverse wave vector of the instability, and V 
corresponds to the relative strength of the self-focusing versus defocusing. Figure |4.1 
shows this for the case of positive V. For our UV problem 

Q/\n 3 

V = fc° 0|£, 2 |£o | 2 = k 0 7i 2 \£ 0 \ 2 + ^onsN", (4.2.13) 

so if n -2 + 1.5 n, 3 1 £q\ > 0, there is an instability for certain values of k±. For the ri 2 
and n 3 that describe UV propagation in air, this cutoff value corresponds to a peak 
intensity of 2.4 x 10 16 W/m 2 . The form of this result is similar to that of Cerda [45] . 
but with our nonlinear response containing both n 2 and terms. The combination 
of the two terms allows the possibility that the system is stable for f! = 0 , unlike the 
result for n 2 only. I 11 the next section, we will examine the stability for nonzero 
and we will obtain the above results as a limiting case. 
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Figure 4.1: Growth rate Re A for spatial perturbations. The growth rate 
is plotted as a function of the spatial frequency parameter T, for the case 
of positive V. For spatial frequencies that correspond to T < 2V, the 
plane wave solution is unstable. 


4.3 Linear Stability Analysis 

To model the stability of long filaments with respect to fluctuations in time and 
space, we abandon the approximation that the plasma density depends only on a 
power of the electric field, and return to using a time dependent evolution equation 
for the plasma. Consider the propagation equation, neglecting losses, combined with 
a simple model for the plasma: 

^ + ik o^\S\ 2 S - (4-3.1) 


dp 

di 


b\S\ 2K — ap 2 . 


(4.3.2) 
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Assume that a steady state solution can be found with £ 0 (r, z) = u(r)e lkz . The steady 
state value for the plasma density will be 


Pair) = 



u(r)\ K . 


Let us investigate a time-perturbed solution consisting of 


(4.3.3) 


£(r, z,t) — £ 0 (r, z) + £+(r, z)e lQt + £_(r, z)e lQ \ (4.3.4) 

P(r, z, t ) = p 0 {r) + p+(r, z)e~ int + p* + (r, z)e m , (4.3.5) 


where £ + , £-, and p + are small compared to £ and p respectively. Since the plasma 
density p is real, we must use p + and p* + as the perturbation coefficients for the 
plasma. 

First, we solve the plasma equation to find p + (r,z), by inserting the perturbed 
solution for p, and keeping only terms linear in the perturbations. We find that 

^ = -ittp + e~ l - lt + i£lp* + e m 
at 

= b\£ 0 \ 2k + bK\£ 0 \ 2K ~ 2 [(£*£+ + £ 0 £*_) e~ m + c.c.] (4-3.6) 

- apo - 2ap 0 (, p + e~ lClt + c.c.) . 

The zero order terms cancel, by the definition of po- Equating terms that oscillate as 

_^'O/ - • 

e gives 


i£Lp + = —2ap 0 p+ + Kb\£ 0 \ 2K ~ 2 (£*£ + + £ 0 £1) , 


(4.3.7) 
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and we find 


P+ = 


Kb\£ 0 \ 2K ~ 2 (£££+ +£ 0 £*_) 

2ap 0 — ifl 


(4.3.8) 


With the terms oscillating as e , we obtain the complex conjugate of this equation. 
Next, we take this solution for p + and the full form for £, insert them into the 


propagation equation (4.3.1), and expand terms, keeping only linear perturbation 


terms. For the zero order terms we obtain 


d£n .o .aur 

~7jp- = V±£o + 'i'koU2\£o\ £o — i 2 Po^tn 


(4.3.9) 


where we assume by construction that £q and po satisfy this equation. For the terms 
oscillating as e~ lClt we obtain 


d£+ 

dz 


—Vi^ + + ik 0 n 2 [2|£ 0 | 2 £+ + £q£1\ 


— i 


ctujt / „ Kb\£ 0 


2K-2 


Po£+ + 


2ap 0 — iQ 


[\£ 0 \ 2 £ + + £ 2 £*_]), (4.3.10) 


and for the terms oscillating as e lflt we obtain 


OS- 

dz 


—+ ik 0 n 2 [2\£ q \ 2 £. + £ 2 £* + \ 


gut ( Kb\£o 

-1—^~ I po£~ + 


2K-2 


2apo + iQ 


[\£o\ 2 £-+ £*£*+]). (4.3.11) 


These two equations give the evolution of the coupled perturbation fields £ + and £_. 
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4.4 Plane Wave Perturbations 

Let us first solve these equations by using an ansatz for the form of £ + and £_ 
that separates the dependence in z from the dependence in the transverse spatial 
coordinate. Recall that So(r,z) has the form u(r)e lkz . If we use the ansatz 

£+(r, z) = u+(r)e Xz £-(r, z) — u_{r)e x * z , (4-4.1) 


the exp(2 ikz) dependence of the £g terms will not be cancelled. However, if we 
modify the ansatz so that the z dependence of £q£1 matches that of £ +1 we will get 
cancellation. To do this, we first make the simplifying assumption that £q, u +1 and 
U- are independent of transverse coordinate. We thus treat the case for which the 
steady state solution is a plane wave, as in the previous section. Let the perturbation 
fields also be plane waves, with the same propagation constant as the £ 0 solution: 


£+(r, z) = u + exp(A z + ikz + ikj_ • r) 
£-(r, z) = U- exp(A*z + ikz — ikj_ • r), 


(4.4.2) 


where k^ is the transverse wave vector k x x + k y y and r is the transverse coordinate 
xx + yy. When we insert these forms into the above equations for £ + and £- and 
cancel the exponential dependence, we find 


(4.4.3) 


(ik + ik _l + A)u + = iB [2u + + u*J\ — iC [poU + + H(P) [u + + u*_] j, 
(ik + ik _l + A *)u- = iB \2u- + u* + ] — iC ^ pou _ + H*(P) [u- + j 


(4.4.4) 
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where 

A(£l) = B = k 0 n 2 \£ 0 \ 2 C=^, (4.4.5) 

2apo — i\l 2 

and the quantity k±, from the contribution from the term, is defined by 


i 

2k 


Vl£ + 




(4.4.6) 


Equation (4.4.3), combined with the complex conjugate of Equation (4.4.4), gives 


a matrix equation 


M 

u + 

= 

ik + ik± + A 


u+ 


u*_ 


—ik — ik l + A 


u*_ 


(4.4.7) 


where 


2 B — (7(/3o + A) 

B-AC 

= i 

X Y 

AC -B 

C\pq + A) — 2 B 


-Y -X 


(4.4.8) 


A nontrivial solution to this is only possible if 


A 2 = Y 2 - (X - k - k ± ) 2 . 


(4.4.9) 


Since k and k\ are both positive, they cannot cancel. As long as A 7 ^ 0, there will be 
both positive and negative solutions for Re A. The above result for A 2 factors as 


A 2 = —{B -Cp 0 -k- k ± )(3B - C{p 0 + 2 A) - k - k ± ). 


(4.4.10) 
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Recall that in the previous section, the propagation constant k for the plane wave 
solution was ko An = |^o | 2 + fco'^ 3 £q | 3 ■ Observe from the definitions of B and C, 

combined with the steady state result for the value of n 3 , that this is also equal to 
B — Cpo for a plane wave solution. Thus we can simplify the result, and present the 
equation for the eigenvalues in terms of and Q. Now we can rewrite 


A 2 = Jc ± [2(B-AC) -k±] 

= k ± (2Y-k ± ) (4.4.11) 

= 2 Y{n)k ± - k\. 


Note that the quantity T from Section 4.2 is equal to k ±, and Y'(0) is equal to V. 


The form of the growth rate is identical to that in Section 4.2 but with V, which 
was a real number, replaced by Y, which is in general a complex function of 12. This 
implies positive and negative roots, and thus unstable growth, for all frequencies 12, 
except where A = 0. 


Figure 4.2 shows an example of the growth rate Re( A) plotted against and 
12. Note that the growth rate is everywhere positive, except for k± — 0 or 12 = 0, 
since we must choose the positive root of A 2 . We expect zero growth for k± = 0, 
since this corresponds to a uniform plane wave intensity correction. Additionally, we 
have chosen an intensity large enough such that the plane wave solution is stable for 
= 0. The form of this growth rate is similar to that obtained by Bian et al. |57j for 
a comparable propagation equation. 

To investigate the behavior for large values of fcj_, let Y = Y' + iY", and make the 
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Figure 4.2: Growth rate vs. k± and fl, for intensity 3.2 x 10 16 W/m 2 . 


following approximation: 


A 2 = k±(2Y — k±) = k 


2 , 2 Y' + 2 %Y" 


A = ±ik i ( 1 


2Y' 2/Y' 


k i 


I 


1/2 


V' 

±ifc, ( 1-^- 


fc_L k_[ 
Re( A) « ±Y" = ±/m(Y). 


(4.4.12) 


From this, observe that the value of Re( A) approaches Im(Y) as k _|_ —> cx). For larger 
values of fl, there will be a corresponding to peak growth rate. However, for small 
values of the peak growth rate is at k —» oo. The value of k± is effectively limited 
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by its relationship to k±, the transverse wave vector, which is assumed to be small 
compared to k z in the paraxial approximation. For a grid with Ax = 0.8 /im, the 
maximum value of k± is approximately 3000 m 1 . 


4.5 Response for Infinite Temporal Frequency 


As —> oo, the instability growth rate does not die off, nor is there an upper 
cutoff, as in some cases for modulational instability that only consider self-focusing 
and GVD [351 53J. As —> oo, Y (Q) —> B , so A 2 —» k±(2B — k±). This limit is shown 


in Figure A3 For values of k± within the unstable range 0 < k± < 2 B, the solutions 
are unstable even as fl —> oo. For the 3.2 x 10 16 W/m 2 solution used to generate the 
previous figure for the eigenvalues, the value of 2 B is 126 m -1 . 

This result is at first puzzling, but it can be explained by examining what happens 
to the equations in the stability analysis. If we take > oo, we remove all the effects 
of the stabilizing 713 or plasma term from the coupled equations, and we are left with 
the stability of the Kerr medium, which is known to be spatially unstable [HU 341 14H1J . 

We also expect difficulty for any numerical analysis of this instability, since the 
frequency response to the instability is unbounded. This means that all frequencies 
present in the numerical model will experience growth, and the growing held will 
not be bandwidth limited, as is desired for a spectral code. However, we will see in 
Chapters [5] and [ 6 ] that we can still observe the initial behavior of the instability by 
using a bandwidth-limited noise seed. 
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Figure 4.3: Growth rate in the limit of > oo, for the £ 0 solution with 
peak intensity 3.2 x 10 16 W/m 2 . 


4.6 Analyzing the Stability of UV Filaments 


The analysis so far suggests three ways to approach the stability problem for UV 
filaments. We have done a plane wave stability analysis, which predicts the growth 
rate as a function of a single transverse wave vector ki and temporal perturbation 
frequency fl. This analysis predicts unstable growth for all frequencies fi. 

Next, we can solve the equations for the coupled perturbation fields in the linear 


stability analysis, given by Equations (4.3.10) and (4.3.11). Section 5.3 describes the 


numerical propagation approach to solve these coupled equations. In this case, we will 
obtain results for all transverse wave vectors that are represented by the numerical 
grid, instead of a single transverse wave vector. However, we are still limited to a 


single temporal frequency fl. The coupled equations can be solved for an arbitrary 
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steady state solution £ 0 . If we use a plane wave solution for £ 0 , we obtain a way to 
verify this propagation scheme against the plane wave analysis. 

Finally, we can model the full three dimensional propagation, by including the 
time dimension of a pulse as the third transverse dimension in the split-operator 


propagation, which will be described in Section |5.4[ This will allow observation of 
unstable growth for all frequencies hi and that we can numerically represent. 

Chapter [5] presents details of the numerical approach taken to propagating the 
wave equation for the single and uncoupled 2D field case, as well as the full 3D case. 
Chapter [6] presents the results of these computations. For consistency, throughout 
the stability results, we will use the steady state solution £ 0 that corresponds to a 
peak intensity of 3.2 x 10 16 W/m 2 . This solution has roughly 470 MW of power. 

We choose this solution because, as discussed above, in the plane wave analysis, 


its intensity is above the cutoff described in Section 42, and thus it is theoretically 
stable for all at D = 0, unlike the 200 MW solution. We will find, however, that 
the general conclusion—that the long pulses suffer from modulational instability, does 
not depend on the specific power chosen. 
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CHAPTER 5 


NUMERICAL MODEL 


5.1 Introduction 


In Chapter [3j we discussed 2D propagation results that verify the form of the 
predicted oscillatory beam solutions for the case of time-independent beam propaga¬ 
tion. To study the stability results from Chapter [4j we need to run the coupled 2D 


and full 3D propagations described in Section |4.6[ This chapter discusses how these 
numerical propagations are prepared and performed in each case, and some of the 
pertinent numerical issues. 

The propagation code is written in Fortran 90, and uses the OpenGL graphics 
library to provide real time 3D graphics feedback as the propagation runs. There are 
three different choices for the form of the held to be propagated. All three choices 
use the split-operator method of propagation, which is detailed in Appendix [Bj The 
first choice uses a single 2D spatial held, on a grid with two space dimensions, which 
we call an XY grid. This models a CW beam with an arbitrary spatial prohle, and is 
used to produce the results in Chapter |3j The second option propagates two coupled 
2D helds with the same XY grid, and is used to evolve the two perturbation helds in 
the stability analysis. The third option defines a 3D held over a grid with two space 
dimensions and one time dimension (labeled XYT). Full 3D simulations are run using 
either the case of periodic boundary conditions in time, or a pulse of finite length. 
The temporal periodic boundary condition case is intended to model a section of a 
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long pulse or CW beam. In this scenario, we are interested in the initial growth of the 
modulational instability over a section of the pulse. In the finite duration pulse case, 
we examine the combined effects of everything in the model, including the instability, 
losses, and pulse edge effects. 


5.2 Propagating a Single 2D Field 

The single 2D held code solves the problem of CW propagation of the equation 


d£_ 

dz 


—V\£ + ik 0 n 2 \£\ 2 £ + ik 0 n 3 \£\ 3 £ + ik 0 ri 4 \£\ 4 £ 


2k 


pW 


\£\ 


2K-2 


£ 


a 


-(l + iur)p£. (5.2.1) 


The n 3 and 714 terms are present for testing and generality. By using the n 3 term the 
plasma density p can be set to zero, and that term ignored, resulting in a simpler 
calculation, if we are working under the approximation that the plasma acts as an 
effective n 3 term. Alternatively, if we wish to use a different plasma model we can 
use the plasma term and set n 3 to zero. The plasma evolution equation is 

^ = b \£\ 2R ~ap 2 + cp\£\ 2 , (5.2.2) 

where b = (3^ / Khu. The third term, which describes avalanche ionization, did not 
significantly affect the results when included. This agrees with the Schwarz-Diels 
prediction that its importance is small for sufficiently short pulses. 

The initial condition for the single 2D held case is a Gaussian beam. This held 
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£i n is defined on a rectangular grid. The user specifies the grid spacings Ax and Ay 
and the number of points in each direction N x and N y . In practice, we always use 
a square grid (N x = N y , Ax = Ay). The initial beam is characterized by its width, 
and either the maximum intensity or total beam power. The equation for the beam 

is 

E in (x, y) = S 0 exp (—(a; 2 + y 2 )/wl). (5.2.3) 

The beam is initially defined with no curvature, so it is at a waist. From the definition 
in Schwarz and Diels ini. for a Gaussian beam profile I = / 0 exp(—2 r 2 /w 2 ), the beam 
power is P = nw 2 Io/2. Thus if we choose to characterize the beam by its power, the 
initial value of S 0 is chosen by setting 

/ 2 P 

£o — = \ —2 • (5.2.4) 

V ^ w o 

By defining £ 0 in this manner, we use intensity units, which is consistent with the 
treatment in Chapter [2] and Appendix [Aj The code stores the field amplitude £ as 
a quantity relative to the initial peak value of Sq. The initial peak intensity is then 
used as a normalization factor to compute the held intensity \S\ 2 when needed. 

The details of the split-operator method, which is used to propagate the held, can 
be found in Appendix |B} In addition, we use absorbing boundary conditions (ABC) 
for the edges of the numerical grid. This helps prevent spurious contributions from 
held energy “wrapping around” the edges of the grid, due to the periodic nature of 
the Fourier transform. The absorbing boundary conditions take the form of a loss 
term in the potential, described by a Gaussian centered at the edges of the grid, so 
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that the absorption is smooth. However, the conditions are not perfect, and there is 
always a small amount of energy that reflects from them, which can crop up in long 
distance propagations. 


5.3 Evolving Coupled Perturbation Fields 


For the case of the linear perturbation analysis, we need to solve the coupled 


equations from Section 4.3 


= ^T V i^+ + ikorki [2|£ 0 | 2 £+ + £%£!] 

.ctujt f Kb\£< 

i —~— I Po£+ + 


0 


1 2K—2 


2apo — id 


[\£ 0 \ 2 £ + + £ 2 0 £*_\y (5.3.1) 


d£*_ 

dz 


-^rvis: - ik 0 n 2 [2|£bi 2 e: + f 0 2 f+] 

Kb 


a LOT 


Po£l + 


2apo — id 


[\£o\ 2 E*_+ £*£+} . (5.3.2) 


for £ + (x,y,z) and £*(x,y,z). These two fields represent perturbation amplitudes 
for a possible modulational instability of the “steady state” held £ 0 , at temporal 
frequency P. The plane wave analysis gives an idea of the growth rate as a function 
of P and the transverse wave vector k±, but with these equations we hope to observe 
the spatial dependence of the perturbation, for one particular frequency P. 

The same basic split-operator approach from Appendix [B] is used for the coupled 
equations. We will define T and V to be matrix operators obtained from the coupled 
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equations. However, we make one simplification first. Note that the value £q is 
present in the growth equations for the perturbation fields. We know that So has 
the form u{r) exp (ikz), where k is the eigenvalue of the steady state problem. In the 
terms with |£o| 2 , the phase factor is no longer present, but it still appears in the Sq 
terms. As we did in the linear perturbation analysis, we modify the form of the fields 
to include the k dependence: Let 

£+(r, z ) = u + (r, z)e^ z £_(r, z) = w_(r, z)e l * z . (5.3.3) 


When we substitute the new fields into the growth equations, the phase from £q is 
cancelled by the phase from £*_, and we cancel an overall phase exp (ikz). We are only 
left with terms involving | <£7 0 1 2 , which will simplify the numerical evaluation. There is 
also a new term iku +1 which comes from the z derivative. 

-iku + H—-V^u_|_ + iB [2 u + + u*_] — iC ( poU+ + A [u + + ul]) , (5.3.4) 

ZirC 

iku*_ — —V]_u*_ — iB \u + + 2 u*_\ + iC (po M _ + A \u + + ul]) . (5.3.5) 

ZirC 


du + 

dz 

du*_ 

dz 


where A, B , and C are dehned as in Equation (4.4.8). Thus the evolution of 


U = 


u+ 


u 


* 


(5.3.6) 


is given by 


dU 

dz 


iTU + iVU, 


(5.3.7) 
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where 


T = 


Vi/2 k 


—Vi/2 k 


(5.3.8) 


and 


V = 


2 B — C (p 0 + A) — k 

B- AC 


X Y 

—(B - AC) 

-(2 B - C( Po - A) - k) 


—Y -X 


(5.3.9) 


Here we have modified the definition of X to be 2 B — C(po + A) — k. In the code, the 
fields u + and w_ are stored in a 3D array with size 2 along the third dimension. To 
apply the exp(iT Az/2) term, we transform to the spectral representation, multiply 
the u + component by the precomputed linear propagator array, and multiply by 
the complex conjugate of the linear propagator array. The linear propagator array is 
the same as that for the case of the single 2D field in Appendix [Bj This is possible 
because the T term is diagonal in the spectral representation of the {u + ,u*_} space 
that we defined. 

However, the exp(iV Az) term is not diagonal with respect to {u + ,u*_}. The 
matrix V is not Hermitian, so we cannot use a unitary basis transformation approach. 
Instead, we directly compute the matrix exponential, using the spatial representation 
for the fields. This is done for every point in the (x, y) grid; note that X and Y 
depend on the coordinates x and y through £ 0 . Additionally, X depends on fl, the 
temporal frequency of the perturbation. Fortunately, we only need to compute this 
matrix exponential once, before the propagation begins, since it does not change; this 
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is another advantage of factoring the k dependence out of the perturbation fields. 

To compute the matrix exponentials, we use a result derived with Mathematica. 
The special form of the matrix makes the result simpler than the general case of a 
2x2 matrix exponential. If 


a 


m = 


-b 



(5.3.10) 


then 


exp(m) 


2s 


(i a + s)e s 


- — (e s 
2 A 



2s 


— (e s 

2s 


(a + s)e 



(5.3.11) 


Here s = \Jd 1 — b 2 . To obtain the result we need, we let a = iXAz, b = iYAz. We 
then store this 2x2 matrix exponential at every spatial grid point. During the split- 
operator propagation, the action of exp(WA^) is obtained by matrix multiplication 
for every point on the spatial grid. 

To perform the propagation, we need to know the value of the the propagation 
constant k for a given solution £ 0 . We can use the definition of the steady state field 
to get a numerical value for k. We know that 


d£o 

dz 


ik£ 0 = Q£ 0 , 


(5.3.12) 


where Q is the evolution operator Q = iV]_/2k + ikori 2 \£\ 2 H-. Multiply both sides 
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of the equation by £q and integrate to find 


ik / |^o| 2 = (QS 0 )S i 


o • 


(5.3.13) 


We can evaluate both integrals numerically, and thus find the value of k from the held 
itself. This provides a numerical check for the value of k, since it can be compared to 


the value that was used to numerically generate the solution for £ 0 in Section 3.4 


The initial condition for the 2D coupled case is to set £ + and £*_ to fields of random 
noise. The spatial noise held is defined as 


£+ = N 0 {l + RA), (5.3.14) 

where A is a uniform complex random value with real and imaginary parts between 
+1 and —1, and N 0 and R are amplitudes. The same process is used to generate 
£-. We apply a hyperbolic tangent window to the spatial noise held, then Fourier 
transform the noise to apply a similar spectral window. The spectral window should 
be relatively smooth, so that sharp edges do not contribute to any sharp structures 
in the held, thereby inducing numerical instability. 

When the program is run, we watch for a shape to grow out of the noise. The 
program keeps track of the maximum amplitude of the perturbation held, and plots 
it on a log scale plot. Once the shape has settled down, this log plot approaches a 
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line, and the code estimates the slope of the line, which is the growth rate Re( A): 


|£+| oc exp(i?e(A)z) 
log \£ + \ = const + Re(X)z. 


(5.3.15) 

(5.3.16) 


To test the validity of the coupled perturbation code, we try solving the coupled 
equations for the case of a plane wave £q as well as a plane wave seed. Using a plane 
wave £q means that |£ 0 | 2 is effectively a constant over the held. The plane wave 
seed confines the initial perturbation to one value of k±. The growth rate obtained 
numerically in the code matched the growth rate i?e[A(f2, k±)} predicted by the plane 


wave analysis in Section 4.3 to within a few percent, for several test values of Q and 
k±. We included the case of f2 = 0 for intensities both above and below the stability 


cutoff in Section |4.2| This confirmed that the coupled held code correctly modeled 
the plane wave limit. 


5.4 Propagating the Full 3D Field 

For the case of the 3D XYT held, we need to solve the equations 

(5.4.1) 

(5.4.2) 


T = ^j-V 2 ±£ + ik 0 n 2 \£\ 2 £ - T— \£\ 2h 2 £ ~ ^(1 + iur)p£. 
| = b\£\ 2K - ap 2 , 


where £ = £(x,y,z,t ) and p = p(x,y, z,t). We proceed in the same manner as for 
the 2D XY case, except that now we need to calculate the plasma density by solving 
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the differential equation, instead of using an effective 713 or the quadratic solution for 
the steady state plasma. 

We use the split-operator approach with 



V = k 0 n 2 \S\ 2 - ^~\S\ 2K 2 - |(1 + iujr)p. (5.4.4) 


The field £ is Fourier transformed to spectral space, where half of the T operator is 
applied, then Fourier transformed back to real space, where the V nonlinear index 
change and loss terms, as well as any spatial absorbing boundary conditions, are 
applied. Finally, we transform to spectral space for the second half of the linear 
propagator. The T operator contains no time dependence, so it is independent of 
time slice in the full grid, and it is applied in parallel to the different time slices. 

To compute the plasma density at a given z position, we numerically integrate 
the plasma equation. One approach is to use a modified Euler method: 


Pn +1 


Pn +1 


p n + Atb\£\ 2K 

Pn +1 


1 + aAtp' n+1 ' 


(5.4.5) 

(5.4.6) 


However, this sequence does not converge to the steady state value 


p = \AA>|s| ,f , 


(5.4.7) 
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unless the product ab \£\ 2h At 2 « 1. This effect was observed when we started with 


the numerically constructed So solution from Section T4 The plasma generated by 
this method was roughly one-third of what it should have been, for At = 3 x lCff 11 
seconds. To obtain an accurate value for the generated plasma within a few percent 
requires At = 3 x 10“ 13 seconds. Alternatively, one could try the sequence 


Pn.+l ~ 


p n + At b\S\ 2K 
1 + aAtp n 


(5.4.8) 


This converges to the correct value. This may also be generalized to include the 
case of avalanche ionization, though in the computations presented here, we neglect 
avalanche ionization. It remains to be tested whether larger time steps can be used 
with this method. One might also consider higher order, but explicit, numerical 
methods. 

The code uses a third method, based on the exact solution to the differential 
equation describing the plasma. This exact solution applies to the case where the 
source term is constant in time, that is, \S\ 2 is constant in time. The exact solution is 
then combined with a piecewise constant approximation for the field. For each step 
let the source term be constant over the duration of the step. Here we will change 
notation, and temporarily redefine b = b\£\ 2K for simplicity. The equation 


= b — ap 2 


dp 

dt 


(5.4.9) 
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has the exact solution 


p(t) 


Y - tanh[v / a6(i — to)], 


(5.4.10) 


where to is a constant of integration. Define p 0 = y/bja; this is the steady state 
plasma density. If we impose the condition that p(t n ) = p n , we find that 


p(t n + 1 ) = Po tanh 



t n ) + tanh 1 



(5.4.11) 


This scheme presents difficulty if p n > p 0 for some n during the iterated solution, 
since the tanh -1 function then becomes complex. However, we can use an identity to 
expand this expression: 


tanh(tanh 1 x + tanh 1 y) = ———. (5.4.12) 

1 + xy 

Then define At = t n+ 1 — t n . We have 


p(tn+ i) = Po tanh 


tanh 1 ^tanh \fab At'j + tanh 1 


— Pn+ 1 — PO 


tanh y\fab Atj + p n /po 


1 + tanh (\fab Atj p n j p 0 


(5.4.13) 

(5.4.14) 


This approach removes the need to use any inverse hyperbolic tangents, and allows 
the perturbed p n to be larger than p 0 . 

Given the above scheme to compute the plasma density, we then solve it for 
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each value of x and y over time. For this reason this portion of the code is easily 
made parallel, since there is no interaction, such as diffusion, between the plasma for 
different values of x and y. If the pulse is confined in time, so that the field is zero at 
the boundaries, then we solve the equation from the low boundary of the time grid 
to the high boundary. However, if we are running the temporal periodic boundary 
conditions case, we iterate the solution, wrapping back around the time boundary 
several times, to ensure that the plasma is continuous over the boundary. In practice, 
three or four iterations are sufficient to verify this. 

For the case of periodic boundary conditions in time, the initial condition is con¬ 
structed by first loading in the numerically generated £$ spatial profile described 
earlier. Each time slice in the full XYT array is initialized with this profile. Random 
multiplicative amplitude noise is then added to this full field. Finally the noise is 
filtered in the frequency domain, both spatially and temporally. The upper limit 
for the temporal noise frequency is given as an input parameter. It is important to 
ensure that this upper limit is an absolute value and does not depend on the grid 
spacing, or frequency grid size. Otherwise, the observed growth will depend on these 
grid parameters. This is clue to the presence of unstable growth for all frequencies Q. 

5.5 Numerical Considerations and Diagnostics 

For the case of XY propagation, one quantity that is measured for each step is 
an estimate of the beam width. This can be done with a second moment technique 
as follows: Take a ID slice of the beam through the origin. Define the width of the 
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Gaussian by w in the expression g(x) = exp (—x 2 /w 2 ). For this form of g(x ), 


fZ* 2 9(x) 2 W 2 
I-oo ff ( X ) 2 4 


w = 2\ 


' J x 2 g(x) 2 
J g(x) 2 


(5.5.1) 

(5.5.2) 


Both of these integrals are easily performed as numerical summations, and the width 
estimate w is updated for each step of the propagation. The code also updates the 
user with the peak held intensity, and an update of the what the width would be for 
a Gaussian beam propagation under diffraction only. 

Diagnostic tests for the 3D propagation case included initializing the spatial held 
with the numerically constructed eigensolution of the propagation, and starting with 
a perfectly hat temporal prohle. For this we imposed periodic boundary conditions 
in time. When this is run, we observe that its amplitude remains stationary, at least 
for the initial few meters, before tiny numerical errors begin to seed unstable growth. 

Another general test for the 3D code involves using only ri 2 and n-j nonlinearities 
rather than modeling the plasma. This results in the beam or pulse being uncoupled 
in time, since all the time slices act independently, but gives a test that can be 
compared to the 2D code. 

The grid spacings are chosen to minimize computational effort while best repre¬ 
senting the numerical helds. During the testing and generation of results for 2D and 
3D propagation, the results were verihed for different values of Ax, At, and Az to 
check the convergence. Section B.2 in Appendix [B] discusses the relationship between 


the grid spacings and propagation step when using a split-operator method. Table 5.1 
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provides the parameters for the numerical grids used in most of the simulations in 
this work. 

Table 5.1: Typical numerical grid parameters. Those used to generate the 
results in this dissertation are given in bold. 


Parameter 

Name 

Value 

Units 

Number of X or Y points 

N X ,Ny 

128 or 256 

— 

Number of T points 

N t 

128 - 1024 

— 

X or Y spacing 

Ax , Ay 

5.0, 8.0 ,10.0 

/irn 

Z step 

A z 

0.5, 1.0 

mm 

Time spacing 

At 

0.5, 1, 2 

IQ’ 12 s 
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CHAPTER 6 


NUMERICAL RESULTS 


6.1 Coupled Perturbation Fields 


This chapter presents the numerical results of the latter two methods outlined in 


Section |4.6| the coupled field propagation, and the full 3D propagations. We first 
examine the results for the propagation of the coupled perturbation fields £ + and S- 


described in the linear analysis in Section [4~3| The coupled field propagation requires 
less computational effort than full 3D propagation, since it involves 2D fields. The 
results tell us about the shape and growth rate of the perturbation fields. However, 
the trade-off is that these results are linear approximations, and we must choose one 
particular temporal frequency P for which we investigate the instability. We seed 
the perturbation fields with random noise, and propagate them using the method 
described in Section 15.31 

We must choose the steady state held £ 0 as an input for the computation. This 


is done by using the computed exact steady state solution from Section 3.4 


corre¬ 


sponding to the given peak intensity or beam power of interest. Alternatively, we can 
choose £ 0 to be a plane wave with a specific intensity, for a verification of the coupled 
held propagation code, discussed in Chapter [5j In this case, the perturbation held is 
initialized with a single value of . 

For the case of using a steady state beam solution So, the results indicate a 
growing instability for all frequencies P > 0, as predicted. The growth rates do 
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not exactly correspond to the plane wave prediction, because the So beam solution 
contains a combination of plane waves with differing k± values. The intensities of these 
plane wave components arc lower than the peak intensity of the beam solution; lower 
intensities tend to correspond to lower growth rates. Depending on the frequency 
D that is considered, as well as the power of the steady state solution, the growing 
instability fields take different shapes. We discuss the results here for D = 0,10 11 , 


and 10 12 s' 1 . We use the steady state solution So computed in Section T4 with a 
peak intensity of 3.2 x 10 16 W/m 2 , corresponding to a power of roughly 470 MW. 
For D = 0, one might not expect any instability growth, since we chose the peak 


intensity of the steady state solution to be above the stability cutoff in Section 4.2 


However, an unstable mode still grows for 0 = 0, no matter how carefully we choose 
the steady state solution or adjust the grid sizes or spacings. The unstable mode was 
circularly symmetric, with a growth rate A on the order of 10' 2 m' 1 . The explanation 
for this unstable growth is that the solution for S 0 is composed of plane waves that 
can have a lower intensity than the cutoff. Thus a small portion of the held seeds 
growth at = 0. 

Since the original stability analysis is linear, there is no dependence on amplitude, 
and the mode will grow without bound as the coupled fields evolve. This is true for 
all frequencies D; the perturbation fields grow to magnitudes well beyond the validity 
of the small perturbation approximation. However, this approach provides a tool to 
discern the shape of the most unstable mode, as this mode will grow fastest, and 
eventually dominate, as the fields propagate. 

Figure [6J~1 shows the form of the perturbation fields for Si = 10 11 s' 1 , and the exact 
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steady state So field with peak intensity 3.2 x 10 16 W/m 2 , which is our standard 
solution. Figures |6.1(a)| and 6.1(b)| show the absolute value and real part of the 
growing perturbation field S + . The field £_ has a nearly identical shape and is not 
shown. The fields take the form of two lobes; When they initially develop, the two 
lobes may not have the same height, depending on the initial noise seed. However, as 
this strongest mode grows, the two lobes become symmetrical. The early asymmetry 
is due to the presence of other, slower growing modes, from the initial noise seed. 
However, the shape of the fastest growing mode eventually dominates. Since there 
is azimuthal symmetry in the propagation equations, the choice of orientation of the 
two lobes is determined by the initial conditions. The rough textured background 
in these figures is the initial noise seed. These figures show the growing solution at 
z = 0.78 m. 


Figure 6.1(c) shows the combination So + S + + S-. At this distance, S + still has 
a somewhat smaller magnitude than £ 0 - As mentioned, the perturbation will grow 
without bound, but we can perform this example summation while the magnitudes 
are comparable. Note that we have arbitrarily chosen a time origin t — 0 such that 


the full equation for the field plus perturbation, Equation (4.3.4), gives just the sum 


of the fields. As the perturbation evolves, the absolute value approaches a fixed 
shape, but the real and imaginary parts continuously change phase. This makes the 
beam become lopsided, and oscillate in the manner of a “snake” instability. After 
the initial growth of the strongest unstable mode from the noise, the growth rate 
of the perturbation fields is exponential. The growth rate for this nontrivial So at 


H = 10 n s 1 is approximately 13 m , which, as discussed and seen from Figure 4.2 
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is lower than the general neighborhood of linear growth rates, given this value of fi. 
The peak plane wave growth rate for = 10 11 s' 1 is roughly 35 m _1 . 


Figure 6.2 shows the form of the perturbation fields for the case of = 10 12 s^ 1 . 
For this frequency, the growing held initially appeared to have the two-lobe structure 
that was observed for = 10 11 s _1 . However, as the held is propagated, the structure 
changes into the symmetrical form shown in the hgures. Thus for H = 10 12 s -1 , the 
growth rate for the symmetric mode is faster than the rate for the asymmetric mode. 
This growth rate is approximately 40 m -1 , while the peak plane wave growth rate for 
this frequency is roughly 60 m -1 , 

For lower frequencies, the shape of the perturbation appears to be a two-lobed 
structure, while for higher frequencies, the shape of the perturbation begins as a two- 
lobed structure, but then evolves into a circularly symmetric form. In both cases the 
initial form shows some asymmetry, as a combination of modes of varying growth 
rates develops. The asymmetric structure at low leads to a “snake” instability, 
while the symmetric structure at higher H causes a “neck” instability. The cutoff 
value where the mode shape changes from two lobes to azimuthally symmetric is 
not known exactly, but will depend on the peak intensity that defines the zero order 


solution. 
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(a) Absolute value of £ + 



(b) Real part of £ + 



(c) Absolute value of £ 0 + £ + + £_ 


Figure 6.1: Coupled perturbation fields for fl = 10 11 s^ 1 . In (a) and 
(b) the absolute value and real part of £ + are shown. Part (c) shows an 
example of combining £ + + £_ with So, which results in a lopsided beam. 
Here the noise seed is visible as the rough texture of the background in 
the figures. 
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(a) Absolute value of £ + 



(c) Imaginary part of £ + 


Figure 6.2: Coupled perturbation field for 12 = 10 12 s _1 . At this fre¬ 
quency, the fields are circularly symmetric. Here the fields have been 
evolved for several meters and have grown far beyond the noise seed. 
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6.2 Periodic Boundary Conditions in Time 


We now turn to the case in which the time domain is fully resolved, and we use the 
three dimensional propagator. For the first study of full three dimensional results, we 
propagate using periodic boundary conditions in time, using the method described 


in Section A4 This allows investigation of the stability of the central region of a 
hypothetical long pulse, without needing to resolve or consider any transient effects 
based on the pulse edges. We initialize the temporal profile as flat, but possibly with 
some small modulation or random noise. We can initialize the spatial profile to a 
steady state solution or a Gaussian beam. Additionally, we can introduce random 
noise over the full 3D spatiotemporal grid. 

The effect of noise depends on the bandwidth of the noise seed. Recall that the 
stability analysis in Chapter [4] indicates that unstable modes are expected to grow 
for all temporal frequencies. Thus, if the noise seed contains all available frequency 
components, the noise will grow for all frequencies allowable on the grid. If the noise is 
spectrally limited to lower frequency components, the spectrum of the growing noise 
will be approximately cut off at this limit; however, higher frequency components 
eventually appear, due to the coupling between space and time in the plasma response. 

As an example, we initialize the held with the steady state solution for Sq shown in 


Figure 3.6, with peak intensity 3.2 x 10 16 W/m 2 . The noise is set to be random noise 


with 1 percent amplitude, bandwidth limited to 101 < 10 12 s 1 and kj_ < 2 x 10 5 m 1 . 


The results of the propagation for three different distances are shown in Figure 6.3 


The subfigures show plots of the absolute value of the held amplitude taken in cross 





section, for the central (y = 0) two dimensional slice of the full three dimensional 
field. The initial condition and early evolution are not shown, as the amplitude of 


the noise is initially too small to see. Figure 6.4 shows the spectral representation of 


the state reached in Figure 6.3(b) The central peak of the spectral representation, 
which is the zero frequency component, has been cropped in the figure, to show 
the spectrum of the growing perturbations. The spectrum grows strongly within the 
region bounded by |Q| < 2 x 10 12 s -1 . This value was the upper limit for the temporal 
frequencies present in the noise seed. Both symmetric and antisymmetric (neck and 
snake) modulations are observed. 

As mentioned, a problem with simulating the growth of unstable modes on the 
full 3D grid is that we can never represent all the temporal frequencies involved, 


since, as seen in Section 4.5, the growth rate is positive for all frequencies. We can 
only demonstrate the growth for a limited range of frequencies. Additionally, as the 
field evolves, the noise grows until it becomes so strong that it exceeds the numerical 
capabilities of a grid; the time or length scales of the features of the noise peaks can 
approach the grid spacing in space or time. However, this analysis provides a useful 
look at the initial growth of a modulational instability. For the UV filaments studied 
here, it is apparent from these results that the fundamental steady state solution 
suffers from the instability and will fragment in time. Thus, for a pulse with finite 
duration, we expect that the transient leading edge of the pulse will not be able to 
converge to the steady state solution that represents the central flat temporal region 
of the pulse in the Schwarz-Diels model. We model such a pulse in the next section. 
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(c) z = 17.5 cm 

Figure 6.3: Growth of instability for periodic boundary conditions in time. 
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omega [1/s] 


Figure 6.4: Growth of instability for periodic boundary conditions in 
time, in the spectral domain, at z = 15 cm. The central zero frequency 
component has been cropped for clarity. The growing noise stays largely 
confined to the area bounded by the noise seed cutoff. 


6.3 Full 3D Model 

To simulate a finite duration pulse, we must prepare the field with a temporal 
pulse profile that resolves the edges of the pulse. We remove the explicit periodic 
boundary conditions in time, by only using one iteration to compute the plasma 
density. In many of the runs, the initial temporal profile is chosen as a supergaussian, 
defined as follows: 


£{t) oc exp 


[d-F >) 2 


M 



(6.3.1) 


where M is the supergaussian order. This allows examination of a pulse that con¬ 
tains both transient effects at the edges and a central flat region. Such “flat-top” 
temporal profiles have been used in experiments [5H] and numerical simulations [33] 
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to study both the edges and the center of a pulse. However, UV pulses of hundreds 
of picoseconds in length, such as those we consider, have not yet been produced. We 
can initialize the spatial profile to match the exact computed stationary solution, for 
diagnostics, or to a Gaussian beam profile, to represent what might be experimen¬ 
tally generated. We can also set the temporal profile to a Gaussian, by setting the 
supergaussian order M to unity. 


As a first example, consider the evolution of the pulse shown in Figures 6.5 6.6 


The figures show an XT cross section of the full 3D XYT field. It was initialized 


with the exact steady state spatial solution, as constructed in Section 3.4, with peak 
intensity 3.2 x 10 16 W/m 2 . The temporal supergaussian was of order M = 10, with 
t p = 200 ps, giving a pulse of roughly 400 ps in length. As the pulse propagates, it 
develops a series of collapse events on the leading edge. These collapses are seeded by 
the transient shape of the edge, which does not match a steady solution. During each 
collapse, the intense plasma causes the held to dip in amplitude and temporarily form 
a ring, which then seeds the next collapse. The trailing edge independently breaks 
apart into a ring, but does not cause any collapse events. These self-focusing collapse 
events appear to consume the pulse at a roughly linear rate; for the power in the 
example, this rate is roughly 200 ps/m. 

If we initialize the spatial profile of the pulse to a Gaussian profile, which makes 
it differ from the exact spatial solution, the beam sheds a small ring in the central 
region, as observed in the 2D results. The overall structure of the collapses, however, 
is not significantly different from the previous case. Figure |fT7| shows an intermediate 
step in such a propagation, and the ring is visible as a bump along the side of the 
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Figure 6.6: Propagation of a finite duration pulse, continued. Note the 
series of collapse events that consumes the pulse. 
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pulse. The initial condition was a spatial Gaussian with width 120 /im. 



Figure 6.7: Propagation of a spatially Gaussian pulse, at z = 25 cm. 
Note the small ring that is generated, as the field sheds power to move 
toward an equilibrium solution of the 2D propagation equation. 


This consumption of the pulse by collapse, while it arises from the transient be¬ 
havior at the beginning of the pulse, does not appear to be itself a transient effect. 
That is, it does not decay—it maintains a steady rate of consumption of the pulse. 

It is unclear exactly what shape the collapsed peaks take. They tend to collapse 
to dimensions comparable to the grid spacing in time. With no loss, there are only 
numerical dispersion effects induced by the finite spectral bandwidth of the grid that 
act to limit the collapse [33]. The simulation has been done with multiple values of At 
to check the convergence, and in all cases the overall result is consistent with respect 
to the spacing of the pulses and the consumption rate. An example of this comparison 
is shown for the on-axis amplitude at z = 25 cm in Figure |(u8| Losses or other physical 
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effects such as GVD will stop the temporal collapse; however, numerically resolving 
these effects requires shorter (femtosecond) time grid intervals, making it impractical 
to model a long pulse. Our result is not the precise shape of the collapse events, but 
their presence, which confirms the unstable nature of the leading edge of a pulse. 



Figure 6.8: Comparison of collapse events for multiple values of At. The 
plot shows the on-axis amplitude of the field at z = 25 cm during the 
initial portion of the pulse, for three different values of At. 


At the farthest propagation distance shown, in Figure 6.6(c), there is noise growing 
on the trailing edge of the pulse. This is evidence of the growth of the modulational 
instability predicted in Chapter [4] and observed in Section C2 The instability has 
grown from numerical noise to a magnitude large enough to be visible in the figure. 
In the next section, we further examine the instability by seeding the initial pulse 


with controlled noise. 
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6.4 Full 3D Pulses with Noise Seed 


Figure [679] shows the propagation of a pulse with the same initial characteristics of 
that in the previous section, but now with a noise seed added to the initial condition. 
The initial condition is not shown, since the initial noise amplitude is small enough 
that one cannot see it initially. The initial noise amplitude was 0.01 percent, limited to 
|P| < 10 12 s -1 . The edge-induced collapse events proceed identically at first, but after 
a few centimeters the noise growth appears, and around z = 35 cm, it overwhelms 
the pulse. The noise growth destroys the pulse in a global fashion, as opposed to the 
collapse events, which consume the pulse from start to finish. 

One might ask if the impact of the modulational instability depends on the tem¬ 
poral shape of the pulse. Does it require a flat, steady-state region? To answer this, 


Figure 6.10 shows one step in the propagation a pulse with Gaussian initial temporal 
profile. Here the width was t p = 160 ps, and the initial noise amplitude 0.04 percent. 
The spatial profile was initialized to a Gaussian with peak power 500 MW and width 
120 fin l. This example demonstrates that the Gaussian pulse still suffers from front- 
edge collapse, as well as the growth of the modulational instability. The particular 
shape of the input pulse does not appear to have a significant effect on the collapse 


or instability. 
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Figure 6.9: Propagation of a finite duration pulse with a noise seed. The 
beginning of the collapse of the leading edge of the pnlse is overwhelmed 
by the exponential growth of the noise. 
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(b) z = 25 cm 


Figure 6.10: Propagation of a temporally Gaussian pulse. The collapse 
events and noise growth take place in the same fashion as for a flat-top 
pulse. 
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6.5 Complete Pulse with Noise and Loss 


The previous sections considered the propagation of pulses under simplified condi¬ 
tions, to isolate the different effects. Let ns now prepare a pulse with Gaussian initial 
spatial shape, a noise seed included, and absorption losses included in the model. We 
will retain the flat-top supergaussian shape in time. These conditions are intended to 
be a better representation of what one might be able to generate in the laboratory. 

The initial condition was a spatial Gaussian defined with peak power exactly 
500 MW, with width 120 /mr, and a temporal supergaussian of order 10 and width 
t p = 200 ps. This corresponds to roughly a 400 ps pulse. The noise seed was 0.01 


percent, bandwidth limited to |P| < 10 12 s 1 . Figures 6.11-6.12 show the propagation 


of the pulse. Additionally, Figure |6.13| shows the induced plasma density at certain 
propagation distances. The characteristic decay time of the plasma can be seen on 
the trailing edge. The plasma density does not need to decay to zero at the time grid 
boundary, because it is an input to the nonlinear response, and is multiplied by the 
field, which is zero at the boundary. 

In all the simulations presented in this dissertation, we use the same base steady 
state solution with power of 470 MW, or a Gaussian beam cross section with power 
500 MW, both of which give a peak intensity of 3.2 x 10 16 W/m 2 . As discussed earlier, 
these powers were chosen so that this peak intensity is greater than the plane wave 


cutoff for stability in the G = 0 analysis in Section 4.2 The specific power, however, 


is arbitrary. The instability and pulse consumption shown here are examples of the 
effects which occur for any power level that supports a self-trapped beam. 
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(a) z = 0 cm 



(b) z = 20 cm 



(c) z = 30 cm 

Figure 6.13: Generated plasma density for the realistic pulse. Plot scales 
are different clue to creation of sharp peaks. 
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6.6 Fluence 

In experiments it is difficult to resolve the temporal profile of a fast, high intensity 
pulse. Instead, often the fluence is estimated, for example by measuring the profile 
of a laser burn spot [1] or with a CCD camera [ID]. For a given propagation distance 
z, the fluence is: 


F{x,y,z) = / \S(x,y,z,t)\ 2 dt , 


( 6 . 6 . 1 ) 


where | £ | 2 is in units of intensity. The units of fluence are J/m 2 , and it corresponds to 
the energy density delivered by the pulse. By calculating the fluence for a simulated 
pulse, we can obtain a quantity than can be experimentally observed. We can also 
compute the width of the fluence profile, using the second moment method, and use 
it as a measure of how confined the beam remains as it propagates. To determine 
the combined effect of the linear consumption and the exponential growth of the 
instability, we examine the plots of the fluence, and the width of the fluence. 

Figure |6.14 shows the fluence profiles corresponding to the propagation in the 
previous section of the 470 MW peak power pulse propagated with noise and loss. 
The impact of the pulse collapse on the fluence depends on the fraction of the pulse’s 
temporal width that has collapsed. For this example, at roughly 30 cm, the noise has 
begun to dominate, which begins to significantly affect the fluence profile. Up to that 
point, the fluence profile is not significantly affected by the collapsing front edge. 


Figure 6.15 shows the width of the fluence profile as a function of propagation 
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(g) z = 35 cm (h) z = 40 cm 


Figure 6.14: Fluence profiles for the realistic pulse, at several propagation 
distances. 
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distance for the same case. The fluence width initially decreases, since the initial 
condition was chosen to have a spatial width larger than the corresponding stationary 
solution. This is in agreement with the oscillations predicted by the Schwarz-Diels 
model in Chapter [3j Then the pulse stays confined as the width is maintained, until 
about 30 cm. During this phase, the leading edge of the pulse is consumed by transient 
collapse. The effect that this has on the width of the total fluence is dependent on 
the length of the pulse; for a longer initial pulse, a given consumption rate in ps/m 
affects a smaller fraction of the pulse, and has less of an impact on the fluence profile. 
However, once the exponentially growing instability becomes significant, in this case 
around 30 cm, it quickly spreads the beam width. The noise is seeded uniformly 
along the pulse, and thus acts as a global effect, and its exponential growth quickly 
dominates. After 60 cm, the fluence width stops increasing because it is numerically 
confined by the absorbing boundary conditions. 

If one is only interested in the fluence at the target, the claim could be made that 
the pulse propagates successfully up to the point where the noise growth dominates. 
This point is determined by the initial amplitude of the noise and the growth rate 
of the noise, which is is roughly 10-20 cm -1 , corresponding to a growth length of 
5-10 cm. For this example, the fluence width stays confined until about 40 cm, 
somewhat more than the Rayleigh range on the order of 10 cm. A pulse could be 
prepared to propagate farther in a confined fashion if the intensities were lower and 
the initial noise was controlled, but not for distances on the order of kilometers, as 


was desired. 
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Distance (meters) 


Figure 6.15: Width of fluence profile for the realistic pulse. The initial 
decrease in width is characteristic of the oscillations predicted by the 
Schwarz-Diels model. This pulse travels roughly four Rayleigh ranges 
before being destroyed by the modulational instability. 
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CHAPTER 7 

IMPLICATIONS FOR UV PROPAGATION 

7.1 Length Scales 

To best understand the impact of the various effects, it will help to consider the 
dominant length scales, and determine which are the most significant. The fundamen¬ 
tal length scale for a Gaussian beam is the Rayleigh range, defined as nwl/X, where 
wo is the beam waist radius. ‘Long distance propagation’ is relative to this scale; if a 
beam does not diffract over many Rayleigh ranges, then it can be considered to have 
propagated a long distance. For most of the example beams that we have considered, 
the waist size wo is on the order of 100 /irn. Combined with working in the ultraviolet 
at 248 nm, this gives a Rayleigh range on the order of 10 cm. 

The length scale for linear (Rayleigh scattering) loss is given by 

— = -j m = 2 km (7.1.1) 

cp 5 x 10- 4 v ; 


in the ultraviolet at 248 nm |26]. This length is the power 1/e length for linear 
loss. There are also length scales for nonlinear absorption, which depend on the 
field intensity. This can be seen from the results of propagating the 2D beams in 


Section T3, where the power in the beam drops roughly exponentially down to below 
200 MW, but the loss rate then slows, because the absorption is nonlinear. This 


scale therefore depends on the power in the beam. Given the power loss shown in 
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Figure 3.5(b), a rough estimate for the nonlinear absorption 1/e length for power is 
20 meters. 

We have both predicted and observed the length scale for instability growth. The 
growth rate A is an exponential growth rate, and so 1/A gives a length scale for e 1 
growth. Typical values for this growth rate are between 10 and 20 m _1 , making the 
growth length scale roughly 5-10 cm. These growth rates become smaller as beam 
power decreases. 


The transient collapse of the beam observed in Section ^3 takes place on a roughly 
linear scale. For the power in the example, this consumption rate was 200 ps/m. This 
means that, neglecting instabilities, for a pulse on the order of 1 ns, the propagation 


distance before complete consumption and collapse is roughly 5 m. Section 6.6 
amined the impact these effects had on the fluence profiles. 


ex- 


Table 7T summarizes some of the pertinent length scales for ultraviolet propa¬ 
gation. Note that, not counting the large uncertainty in the Raman growth, the 
instability growth occurs over the shortest distances. The instability growth scale is 
comparable to the Rayleigh range. The absorption, both linear and nonlinear, is not 
strong enough to counter this instability before it fragments the pulse. 


7.2 Experimental Implications 

Given the numerical results from Chapter [6j we expect that this instability and 
transient consumption of the pulse will have a significant effect on the propagation 


of ultraviolet filaments. 
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Tabic 7.1: Typical length scales for 470 MW peak power UV filaments 


Name 

Value 

Units 

Rayleigh range 

10 

cm 

Linear loss length 

2000 

m 

Nonlinear absorption length 

ps 20 

m 

Instability growth length 

ps 5 - 10 

cm 

Pulse consumption distance 

t p /(200 ps/rn) 

m 

Raman scattering 

0.01-10 

m 

Upper limit (Avalanche ionization) 

~ 4-60 

ns 

Lower limit (Plasma rise time) 

ps 30-200 

ps 


For our example case, a UV filament of length on the order of a nanosecond may 
propagate for distances on the order of a meter before it is consumed by collapse 
events. As discussed in the previous section, this still may be described as long 
distance propagation, as it is many Rayleigh ranges, given the spatial width of the 
filament, ffowever, the results do not support hopes that kilometer scale atmospheric 
propagation of a stable pulse in the ultraviolet will be feasible; the proposal that 
simple UV pulse results can scale to longer durations is not supported. For the 
purposes of energy delivery, one may not be concerned with fragmentation of the pulse 
in the time direction, as long as the fluence profile is confined; these results do not 
rule out the propagation of some complicated dynamic mode over longer distances. 


ffowever, as we have seen in Section |6.6[ the modulational instability also affects 
the fluence profile by spreading it. The rate at which the modulational instability 
fragments the main body of the pulse is faster than the transient consumption. The 
point at which this fragmentation occurs depends on how much noise is initially 


present, but due to the exponential growth rate, it cannot be realistically suppressed 
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for long distances. 

As mentioned in the Introduction, experimental generation of ultraviolet pulses 
of sufficient length to compare with these results has not yet taken place. However, 
if such pulses were to be generated in laboratory conditions and measured, we would 
expect to be able to compare the qualitative behavior of the width of the observed 
fluence profile with our predictions. 


7.3 Theoretical Implications 


A steady state or long pulse model may simplify theoretical predictions, but if 
the physics of the problem admit the effects of modulational instability combined 
with other transient effects, it is dangerous to apply such a model. The presence of 
such modulational instabilities or transient effects, which we have demonstrated for 
the case of UV filament propagation, violate the basic steady state beam assumption 
that is used in the long pulse model. 

We have shown that UV filaments should suffer from these effects, limiting their 
effective propagation range to meters. This will prevent the steady state theory from 
scaling to longer UV pulses as desired H3 fl6j . One might ask if there are other 
regimes for which we may be able to reduce the impact of these problems. However, 
an initial advantage of choosing the UV was that the effect of avalanche ionization 
was greatly reduced, and thus could be ignored, up to certain pulse lengths. Moving 
to longer wavelengths drastically shortens this upper limit H7j. 


Section 2.3 mentioned that there are many applications of variational methods to 
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nonlinear optical propagation problems. The method of Schwarz and Diels H3. when 
modified, was found to give the same results as the variational method of Anderson 
and Bonnedal [2S] , The drawback of a variational method is that it mathematically 
confines the solution to fit the variational ansatz. If this ansatz does not capture 
a mode of a physical instability that is present, the variational method will not be 
able to model the instability. In the case of UV filaments, the current method is 
unfortunately an oversimplification of the dynamics of long pulse propagation. 

In general, caution must be used when describing soliton solutions of reduced 2D 
time-independent versions of full 3D equations. For example, the treatment by Skarka 
et al. of vortex solitons [30] relies on considering the evolution of a single time slice of 
a pulse, which simplifies the propagation to a time-independent form [15] . This allows 
the use of a variational technique to describe the spatial form of the vortices. They 
emphasize the spatial stability of the vortex solutions, but do not consider the full 
temporal dynamics or stability of the pulse. This time-independent method limits the 
applicability of the claim that such solitons will be able to propagate for appreciable 
distances in air, until a full temporal analysis is performed. 

Similarly, the time-integrated approach of Berge et al. to multiple hlamenta- 
tion 05 cannot describe effects such as dynamic replenishment mm and the col¬ 
lapses in both space and time. The complicated dynamic nature of the problem may 
mean that ” long-range” filaments qualitatively describe the results, but the most 
convincing theory must account for the time dependent characteristics. 
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7.4 Raman Scattering 

The application of a long pulse model to realistic ultraviolet filaments raises the 
question of the impact of Raman scattering. The Raman effect is observed when a 
material scatters light to frequencies different from that of the source, through virtual 
energy level transitions [Tj . The new lower shifted frequencies are labeled Stokes lines, 
while the higher shifted frequency components are anti-Stokes lines. Spontaneous 
Raman scattering is usually weak, but stimulated Raman scattering, which can occur 
for intense laser beams, can convert a large fraction of the incident energy to the 
Stokes frequency, in a conical emission pattern [T|. In air, the dominant process is 
stimulated rotational Raman scattering (SRRS) in nitrogen [33], [59], for which the 
particular rotational transition gives a shift of roughly 1.4 x 10 13 s^ 1 [33] . 

The total nonlinear index 77,2 in the Kerr effect for a CW beam includes contri¬ 
butions from bound electrons as well as vibrational and rotational terms. Each of 
these contributions is due to the response of a material mode (electronic, vibrational, 
or rotational) to a driving electric field. This response is proportional to the field 
intensity. For long pulses, the time scale of the material responses is fast, so their 
time dependence averages out, and the total contribution can be included in a single 
effective instantaneous value for n 2 . For shorter pulses, the time scale may be compa¬ 
rable to the rotational response, so the time dependence, or memory, of this portion 
of n 2 must be considered. 

To include SRRS in a propagation model, the n 2 term is split into time dependent 
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and time independent parts [HHIHE]: 


ik 0 n 2 \£\ 2 £ —> ik 0 {l - f)n 2 \£\ 2 £ + ik 0 fn 2 



dt'R(t-t')\£(t')\ 2 


£. 


(7.4.1) 


Here / is the fraction of the nonlinear response that comes from the time-dependent 
component, in this case SRRS. The function R(t) is the impulse response of the 
rotational mode. 

This impulse response is often modeled as a sinusoidal oscillation multiplied by a 
decaying exponential m Eli ng. The decay models the dephasing due to the fact 
that the true response is due to multiple states in a rotational manifold, each with a 
slightly different frequency. However, this simplified model does not take into account 
the recurrences which occur over longer time periods, on the picosecond scale [ 60] . 
For femtosecond pulses, these recurrences might be neglected, by using the simple 
exponential decay. However, for picosecond-nanosecond scale pulses, they may need 
to be included. 

Penano et al. modeled SRRS in air for long laser pulses. Their model considers 
propagation of pulses under the instantaneous (Kerr) term and the delayed (Raman) 
term. They performed a linear stability analysis to find the growth rates for Stokes 
and anti-Stokes fields. Their result is that the Raman gain coefficient for a CW pump 
is (33] 


90 ~ 2T 2 c 


(7.4.2) 


Here tir is equivalent to fn 2 , the rotational contribution to the total value of n 2 in 
the CW limit. ff 0 = uj\ + T|, where a Jr is the Raman transition frequency; u 0 is 
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the pump optical frequency. T 2 is the damping rate of the exponential in the model 
described above, and is labeled the dipole dephasing rate in Penano et al. [33] 

The experimental value of this gain coefficient, for long pulses at 1 /irn wavelength, 
is 2.5 crn/TW [59] . This value is multiplied by the intensity to obtain a growth rate. 
From the dependence on ujq, we expect roughly a factor of four increase when working 
in the UV, or 10 cm/TW. Given the intensities in this dissertation, on the order of 
10 16 W/m 2 , this predicts a growth rate on the order of 10 cm -1 . This corresponds 
to e 1000 growth over one meter, which contradicts experimental results [4], For the 
current picosecond scale UV pulse experiments, the Raman contribution to n 2 was 
not found to be significant; there was little growth, due to the lack of significant 
spectral broadening before the pulses ionized M. 

This basic result is modified by the coupling between Stokes and anti-Stokes waves 
which occurs for small values of the phase mismatch between the two waves. For 
perfect phase matching, there is zero growth; for small values of the mismatch, the 
growth rate is reduced [[33]. The phase mismatch depends on the propagation angle, 
or equivalently the value of k±. For k± = 0, Penano et ah fold the peak growth 
rate to be roughly one-third of the full g 0 value; it decreases to zero for a particular 
angle, and then rises to approach go for large angles. Thus in general, for a beam 
constructed of a continuum of propagation angles, their analysis predicts unstable 
growth on the order of g 0 . 

The central question is whether or not the simplified model of the Raman impulse 
response is useful when considering long pulses. The recurrences may be too impor¬ 
tant to ignore; they affect any characterization of the response by a time scale. The 
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fast oscillations in the SRRS for N 2 take place at ujr ~ 1.4 x 10 13 s^ 1 . This implies 
that temporal grid resolution roughly ten times finer than we have used may be ade¬ 
quate to numerically model the Raman scattering. However, we need to understand 
how to treat the Raman response before performing this calculation. 
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CHAPTER 8 
CONCLUSION 

8.1 Summary 

We have considered the propagation of long filaments in air, primarily in the ultra¬ 
violet. It had been speculated that previous successful demonstrations of ultraviolet 
filaments would scale to longer pulse durations (rang, allowing the propagation of 
long filaments and delivery of higher energies over long distances through the atmo¬ 
sphere. However, the existing theoretical model for long ultraviolet filaments did not 
address their stability, or how they are created im 

To address these issues, we have considered and extended a long pulse model for 
ultraviolet filaments H7|. We have numerically shown that while it is mathemati¬ 
cally valid for the domain of pulse lengths that it claims, the pulse instability and 
transient effects cannot be neglected. We investigated the question of stability using 
three levels of approximation: the plane wave stability results, the evolution of the 
linearized coupled perturbation fields, and the simulation of full 3D pulses using a 
time-dependent plasma model. The results were consistent on all three levels, and 
showed that the modulational instability ultimately results in the fragmentation of 
the pulse. We also demonstrated that the transient shape of the front edge of the 
pulse induces a series of collapse events, which also tend to consume a long pulse. 
Since the long pulse model can no longer be used for a pulse that is fragmented in 
time, a new model will be needed. 
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In the case of ultraviolet filaments, we still predict that propagation over several 
Rayleigh ranges can be achieved, depending on conditions. The transient collapse 
consumes on the order of hundreds of picoseconds of the pulse per meter of propaga¬ 
tion, but keeps the fluence profile relatively confined. The modulational instability, 
however, destroys the entire pulse as it grows, with a gain length on the order of tens 
of centimeters, depending on how it is seeded. The values of these consumption rates 


depend on the model parameters in Table 3.1 which are not accurately known, but 
the above detrimental effects will still be present even for alternate values. The the 
propagation of long stable pulses described by a steady state model over kilometer- 
scale distances does not appear feasible, given the effects that we have considered. 


As discussed in Section 7.3, the long pulse model and other methods that eliminate 
the time dimension have limited prediction power when temporal stability effects are 
not considered. A steady state model is not suitable for the description of a highly 
fragmented pulse. However, there are still open questions as to other effects that may 
be relevant to the propagation of UV filaments. 


8.2 Future Work 

The stability results point to unstable growth for all temporal frequencies, inde¬ 
pendent of the physical parameters such as the wavelength and plasma generation 
rates. We have applied the results to the UV filament regime, but the specific nu¬ 
merical results are only as accurate as the input parameters. The MPI rate needs 
significantly more study before a more certain value can be obtained; this value affects 
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the effective n 3 as well as the time scales of the upper and lower limits in the original 
model HZ|. The wavelength dependence also invites more detailed study; one would 
like to know how if there is a laser wavelength for which the pulse consumption rate, 
as well as the growth rates of modulational and Raman instability, can be reduced, 
but a long pulse model still applied. 

More sophisticated models of the plasma may be incorporated, including such 
effects as plasma diffusion, or variable MPI rates depending on intensity, or pulse 
length. More study on the impact of stimulated Raman scattering for long pulses 
is needed, including how to describe the Raman response for long times. When this 
is understood, the Raman effect could be added to the numerical model, given an 
order of magnitude increase in computation time or power to accommodate a finer 
temporal grid. 

If one wishes to work with long pulses that will fragment as described, is there a 
way to model the collapse events with an approximate theory or suitable variational 
method? The specific structure of the instability modes may also be investigated. 
Is there some way to prepare a pulse so that the transient collapse effects at the 
pulse front are avoided? Alternatively, is there an acceptable level of fragmentation 
of pulses? Could a collection of multiple filaments, interacting in space and time, 
propagate by replenishing the pulse in time after breakup, as has been seen for spatial 
filaments |41j? 

Study could be done on the impact of temporal instabilities to other cases; for 
example, optical vortex solitons [3D]. One may examine directed energy beam appli¬ 
cations, which may not use such small waists and such high intensity fields, following 
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the example of the case of Raman instability [33]. For such atmospheric applications, 
the seeding of modulational instability from scattering or turbulence may be impor¬ 
tant. Finally, a comparison with experiment is needed for the 3D stability results; 
long UV pulses in the hundreds of picoseconds need to be generated in the laboratory. 
Then a comparison between prediction and experiment can be made for the fluence 
profiles, and propagation distances. 
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APPENDIX A 

COMMENT ON UNITS 

The problem of multiple systems of units must be addressed. In this dissertation, 
we use MKS units; a thorough discussion of how quantities in nonlinear optics relate 
in MKS and CGS (Gaussian) units can be found in Boyd (T., 33j. The convention for 
the nonlinear index terms can be bewildering; Marburger [34] lists twelve different 
conventions for the index change term rz 2 1 •S’ | 2 • To avoid unnecessary confusion, this 
Appendix provides some additional clarifications to the units used in this dissertation 
and how they relate to those used in the paper by Schwarz and Diels nzt We use 
many of the parameters from their paper in our numerical simulations. 

In section IV. B of Schwarz and Diels na. the authors convert the evolution equa¬ 
tion for the beam width from using amplitude units to intensity and then power 
units, to write the results conveniently in terms of the critical power. To do this, they 
convert x 4 t° n 3 and finally the plasma absorption coefficient /Jj. The equations 
they obtain are correct, but there is an error in the assertion below equation (10) in 
their paper, which states that n m -i = eocnoh m _i/2 is the relationship between the 
amplitude and intensity refractive indices. The correct relation can be derived from 
the following: The index of refraction of air is expanded in terms of intensity as 

(A.l) 


n — n 0 + nil + n 3 I 15 . 
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Keeping this form in mind, expand the susceptibility in terms of amplitude: 


n 2 = 1 + x 

(A.2a) 

= l + X (1) + X ,S, |f| 2 + X (4) |f| S 

(A.2b) 

= nl + X (3> |£| 2 + X (4 W 

(A.2c) 

-nS( 1 + *T + *TV 

(A.2d) 

V n o n o ) 


Thus, for small nonlinearities, 


y , x (3) ifi 2 , x (4, iy 3 \ 

"~"T + 2 nl + 2nl ) 

(A. 3a) 

x |3) |£| 2 , x (4) |£| 3 

■” 0+ 2n 0 + 2 n 0 

(A.3b) 

= n 0 + n 2 \S\ 2 + n 3 \S\ 3 . 

(A.3c) 


where, by equation (10) of Schwarz and Diels [T7], x ^ — 2 n 0 n( m -i) defines n m -i- 
Next use the relationship between intensity and held amplitude: 


/ 


^oe 0 c, , 2 
2 1 1 


(A.4) 


to write 


n 


n 0 


-n 2 I + 


1.5 


n 3 I 


1.5 


e 0 cn 0 


e 0 cn 0 


(A.5) 
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By comparing this with Equation (A. 1), we find 


e 0 cn 0 _ 

n2 = —i n 2 , 


/eocno 
V 2 



77m 


fe 0 cn 0 \ m / 2 _ 

V 2 ) Um ' 


(A.6) 

(A.7) 

(A.8) 


This differs from the assertion below Equation (10) in Schwarz and Diels [T7] by 
the addition of the exponent on the prefactor, which depends on the order of the 
nonlinearity. 

To simplify computations, one can use either intensity or held units for the quan¬ 
tities S, 7i2,7i3,..., and , as long as the units are consistent. Intensity 

units for S are defined such that \E \ 2 is the intensity. Thus ri 2 \£\ 2 represents the same 
unitless quantity as long as both factors are measured consistently in intensity or held 
units. 

Note that under strict radiometric convention, the quantity with the units W/m 2 
is denoted irradiance, while intensity is reserved for the quantity measured in W/sr. 
However, many authors in this held use intensity for W/m 2 , so we will do the same. 

Schwarz and Diels HZ] use the notation n 2 for the nonlinear self-focusing index 
in intensity units, and 77-2 for the same quantity in held units. For simplicity, in the 
body of this dissertation, we use only intensity units, and use 712 , 713 ,... to denote 
the self-focusing index and higher order indices in these units. 

Next consider the equation for the nonlinear attenuation of the beam. Equa- 
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tion (17) in the Schwarz-Diels paper [T7j gives this as 


dl 

dz 


-3frwN 0 a [3) I 3 + y n 3 I 15 - 2 / 0 


d(w 2 ) 

dz 


(A.9) 


Here / is the on-axis intensity, hu is the energy per photon, No is the density of 
oxygen, cr® is the three-photon absorption cross section, 77-3 is the plasma nonlinear 
index in intensity units, l is the mean free path length of the electrons in the plasma, 
Jo is the initial intensity, and Wq and w are the initial and evolving beam width. 

There are a few corrections that need to made to this equation. In the first term, 
which describes MPI loss, the units of cd 3 ) are listed incorrectly in Table I of their 
paper; the correct units are m 6 s 2 /J 3 . These correct units are used in section II.B of 
the paper. 

The second term describes plasma absorption, but the intensity I should be raised 
to the power 2.5. This makes it agree with equation (16) of the paper, which states 
that the coefficient for plasma absorption scales as I 1 ' 5 . The units of h 3 are also listed 
incorrectly in the table; they are m 3 /W 3 ^ 2 . 

The third term represents the change in intensity due to the change in beam 
width. It is apparently obtained by differentiating I(z) = ( wl/w 2 )I 0 , but there is a 
factor of two present in the equation in the paper that should not be there. Thus the 
correct equation should read 


dl 

dz 


-3 7kWV„<T l3) / 3 + 1* 5 - /„Ai 

l w{z ) 4 


d(w 2 ) 

dz 


(A. 10 ) 
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Next, they convert the equation for the on-axis intensity to an equation for the 
attenuation of the beam power, using the relation that for a Gaussian beam with 
I(r,z ) = I(z) exp(— 2r 2 /w 2 ) the power is P = ttw 2 I 0 /2. The result is 


1 dP 
P~dz 


-p 3 — t P 2 - ft 
■ur 



(A.ll) 


where 


& = 0) 3hwN 0 a^ 



and /? 4 



(A.12) 


are respectively called the three-photon power attenuation coefficient and the non¬ 
linear plasma absorption coefficient m This result can be derived from Equa¬ 


tion (A.10), so the errors in Equation (A.9) do not affect the results in the remainder 
of their paper, other than requiring some corrections to the units in Table I of their 
paper. Table A.l summarizes these corrections to the units. 


Table A.l: Corrections to units in Schwarz-Diels paper 


Description 

Name 

Corrected Units 

Plasma nonlinear index (intensity) 

n 3 

m 3 /W 3 / 2 

Three-photon absorption cross section 

cr( 3 ) 

m 6 s 2 /J 3 

Plasma nonlinear absorption (intensity) 

Pmpi 

m 3 /W 2 

Plasma nonlinear absorption (power) 

(3 a 

m 2 /W 3 / 2 

Three-photon absorption (power) 

(h 

m 3 /W 2 
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APPENDIX B 

SPLIT-OPERATOR METHOD 

B.l Split-Operator Propagation 

For all three grid types that are considered in this dissertation, we use the split- 
operator method to propagate the fields. In the split-operator propagation scheme, we 
define the propagator as a combination of a linear diffraction piece T, and a nonlinear 
index change piece V. Let 

BE 

— =iT£ + iV£, (B.l) 

oz 

where for the example of single 2D field propagation, 


T = 


2 k(] 


V = k 0 n 2 \£\ 2 + k 0 n 3 \£\ 3 + k 0 n 4 \8\ < 


p( R ) 


\£\ 


2K-2 


O 


-(1 + iur)p. 


(B.2) 

(B.3) 


Note that T will be diagonal in the spectral (k) representation, while V is diagonal 
in the spatial (x) representation for £, and V acts as a changing potential. The 
split-operator method represents the incremental step propagator exp[?'(T + V)Az] as 


exp 



(iVAz] 


exp 


TAz 


+ 0(Az 3 ). 


(B.4) 
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For the T piece, we use the spectral representation: 

TAz 

exp(i^—)£(k x , k y ) = exp(i(-k 2 ± /Ak 0 )Az)£(k x , k y ). (B.5) 

Fourier transform £ to the k representation before multiplying it by the exponential 
exp(— ik'j^Az/Ako). This exponential factor can be precomputed for each point on the 
k grid, that is, each value of k±. Thus the application of the two linear diffraction half¬ 
steps is reduced to applying an FFT to get to the spectral representation, multiplying 
the array by the precomputed linear propagator array, and inverse FFT to return to 
the spatial representation. 

The exponential with V, the nonlinear response, is applied in the spatial repre¬ 
sentation: 

exp(iVAz)£(x, y) = eyzp(iV(\£\ 1 )Az)£(x,y). (B. 6 ) 

First we compute the value of V{\£\ 2 ), using the value of £ at the previous z step, 
then compute the exponential factor, and multiply the held in real space by this piece 
of the propagator. The computation of V can be simple if it only involves n 2 and 
71 , 3 , or complicated if a more complete model is used where the plasma density p is 
needed. 

In the code, the held is initially defined in the spatial representation. An initial 
FFT step is performed to transform to the spectral representation before entering 
the main loop. The main loop consists of applying half the linear propagator, trans¬ 
forming back to x space, applying the nonlinear response, transforming to k space, 
and applying the second half of the linear propagator. The two halves of the linear 
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propagation step can be combined in a practical code, but for rigorous accuracy one 
must keep track of where one is in the propagation and make sure for the final step 
that the full propagator is completed. 

The initial field is defined over a space grid of N x by N y points with uniform grid 
spacing Ax and Ay. The space grid is defined so that the point with space coordinates 
(x — 0, y — 0) is stored at location (N x / 2, N y / 2) in memory. For optimal FFT speed 
N x and N y are chosen as powers of 2. This spatial grid defines a spectral (k-space) 
grid where Ak x = 2n/(N x Ax ) and similarly for A k y . The maximum value of k x is 
given roughly by n/ Ax. Recall that the common convention for FFT storage is to 
locate the point (k x = 0, k y = 0) as the first element in the FFT array, and that the 
array wraps around near the Nyquist frequency [61]. 

B.2 Numerical Issues 

When using the split-operator spectral method, one must choose the numerical 
grid sizes wisely. Both the spatial and spectral grids must contain adequate resolution 
and extent to represent the full held. If the held begins to collapse to scales that are 
too small spatially, the spectral representation will how oh the spectral grid boundary. 
The interactive graphical implementation used in our code significantly aids the user 
in avoiding bad choices of grid parameters. 

For split-operator propagation, there are approximate relationships between the 
propagation step size and the grid sizes to ensure accuracy in the method. This can 
be seen by considering the linear di h r act, ion piece of the split-operator formula in 
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Equation (B.4). The individual portions of the split exponential should be small, and 


thus the arguments of the exponential phase factors should be smaller than unity, 
if the approximation is to be good. Except in the case of extreme self-focusing, the 
T term is usually larger than the V term. For our UV propagation example, take 
self-focusing only to get an upper limit; the value of V on axis is 


V = k 0 n 2 \£f ~ 60 m \ 


(B.l) 


The maximum absolute value of T for a given grid is given, in the spectral represen¬ 
tation, by 

• »' 2 ( VA *) 2 (B2) 


k 2 

| T | ri/ max 


2kn 


2kn 


which gives roughly 3000 m 1 for 248 nm and Aa; = 
in the exponential is then 

\T\Az ti 2 Az 


AknAx 2 


/nn. The maximum argument 


(B.3) 


Setting this to be less than unity gives 


^ Ax 2 

7l z 


(B.4) 


as the condition for Az. For Ax = 8 /irn (see Table 5.1), we hnd that A z < 0.66 
mm is the condition, which is just satisfied by our choice of Az = 0.5 mm. If this 
condition is not well satisfied, artificial numerical instabilities can arise. However, if in 
practice the spectral representation of the field does not sample the outer edges of the 
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spectral grid, one can effectively use a larger propagation step. But when studying 
stability, noise on the entire grid becomes important. This was observed in the case 
of propagating the coupled perturbation fields for = 0; if the value for Az was too 
large, unstable artificial numerical growth occurred. A more sophisticated study of 
the split-operator beam propagation method for the case of a fixed potential term V 
was done by Thylen, who also found that the condition on the step size could be 
relaxed in practice [62j . 
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